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CHAPTER  1 


INTRODUCTION 


1 .1  Background 

Recently  much  work  has  been  directed  towards  developing 
constitutive  models  to  represent  the  complex  load-deformation  behavior 
of  soils.  The  models  developed  to  date  have  primarily  been  for  tne 
special  cases  of  dry  and  completely  saturated  soils.  The  use  of  these 
constitutive  models  in  representing  the  behavior  of  partly  saturated 
soils  has  resulted  in  inaccurate  predictions  of  soil  response.  It  is 
the  intent  of  thi3  research  to  formulate  a  constitutive  model 
describing  the  behavior  of  partly  saturated  soils. 

There  are  essentially  two  approaches  which  have  been  used  by  those 
attempting  to  develop  constitutive  laws  for  3oils.  The  first  approach 
is  termed  pnenomenological  modeling.  Phenomenological  models  may  be 
defined  as  those  concerned  with  describing  material  behavior  on  the 
size  scale  of  the  experiment.  For  soils,  thousands  to  millions  of 
soil  grains  and  pores  would  be  included  in  a  model  representation  of 
this  type.  Pnenomenological  metnods  or  tneones  include  emprical 
curve  fitting,  elastic  theories,  elastic-plastic  theories,  and  visco¬ 
elastic  theories.  These  methods  and  continuum  theories  are  concerned 
with  describing  the  overall  observable  behavior  of  the  soil  mass. 

They  are  not  concerned  with  describing  the  actual  mechanisms  causing 
deformation,  which  act  on  the  level  of  tne  grains  and  pores  which 
comprise  tne  soil  mass.  The  second  type  of  approach  is  termed 
micromechanical  modeling.  This  approach  attempts  to  derive 
constitutive  laws  by  considering  the  deformation  mechanisms  acting  on 


a  very  small  but  representative  sample  of  the  material.  For  soils,  a 
micromechanical  model  description  might  include  one  to  hundreds  of 
grains  and  pores  in  the  model  description. 
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The  primary  problem  with  constitutive  models  representing  soils, 
is  a  failure  to  describe  all  aspects  of  their  load-deformation 
behavior.  While  a  model  may  give  reasonable  predictions  under  one  set 
of  input,  it  may  fail  to  predict  the  soil  resporse  under  another  set 
of  input.  With  the  present  knowledge  it  appears  that  a  constitutive 
model  representing  all  aspects  of  soil  behavior  may  not  be  obtainable. 
This  is  due  partly  to  a  lack  of  understanding  of  the  mechanisms 
causing  soil  deformation  and  partly  due  to  the  mathematical 
complexities  one  may  encounter  when  modeling  soils.  A  micromechanical 
approach  to  the  constitutive  modeling  of  soils  may  provide  a  better 
means  to  understand  the  soil  load-deformation  mechanisms. 

1 .2  Scope  of  Work 

The  purpose  of  the  research  studies  contained  in  this  report  is  to 
develop  a  constitutive  model  representing  the  load-deformation 
behavior  of  soils.  The  following  types  of  investigations  are 
contained  in  this  report: 

a)  Review  of  the  available  literature  on  previously  developed 
constitutive  models  describing  the  load-deformation  behavior 
of  soils. 

b)  Development  of  a  constitutive  model  to  represent  the  response 
of  partly  saturated  soils  using  micromechanics,  under 
idealized  conditions. 

1 .3  Method  of  Approach  and  Organization 

The  studies  undertaken  to  achieve  this  objective  are  described  in 


the  subsequent  chapters. 

Chapter  2  contains  a  brief  description  of  experimentally  observed 
load-response  behavior  of  soils  and  review  of  previously  developed 
constitutive  models  describing  this  behavior. 

Chapter  3  contains  the  solution  to  the  problem  of  an  elastic 
sphere  in  contact  with  a  number  of  neighboring  elastic  spheres. 

Chapter  4  contains  the  development  of  the  effective  moduli  to 
describe  the  elastic  behavior  of  a  number  of  spheres  in  contact  while 
surrounded  by  a  liquid  matrix  phase. 

Chapter  5  contains  the  extension  of  the  work  described  in  the 
previous  chapter  to  three-phase  systems. 

Chapter  6  presents  conclusions  and  recommendations. 


CHAPTER  2 


4 


PREVIOUS  WORK 

2 . 1  Load-Deformation  Benavior  of  Soils 

Wnen  a  soil  mass  is  subjected  to  any  arbitrary  set  of  surface  tractions, 
the  result  is  a  volume  and  shear  deformation  of  the  soil  mass.  The  resulting 
displacement  ana  stress  fields  within  the  soil  mass  depend  on  a  number  of 
variables.  Tnese  vanaDles  include  the  type  of  loads  applied  to  the  soil 
mass,  tne  stress  history  of  tne  soil  mass  and  the  chemical,  and  the  physical 
properties  of  the  soil  mass. 

Experimental  observations  of  the  response  of  a  soil  mass  to  various 
applied  loads  nave  provided  a  great  deal  of  information  concerning  the  load- 
deformation  behavior  of  soils.  Tms  information  provided  by  experimental  work 
will  be  briefly  discussed  in  the  next  four  subsections. 

2.1.1  Soil  Behavior  During  Initial  or  First  Loading 

The  term  initial  or  first  loading  will  be  understood  to  refer  to  a  state 
of  stress,  occurring  within  the  soil  mass,  which  the  soil  is  experiencing  for 
the  first  time.  In  Fig.  2.1,  the  initial  loading  curves  correspond  to  those 
lying  beteeen  points  1  and  2,  and  points  3  and  4.  The  arrows  shown  in  Fig. 

2.1  indicate  the  load  path  taken. 

The  amount  of  volume  deformation  resulting  from  an  increase  in  stress  will 
depend  on  the  relative  density  of  the  soil  mass.  The  relative  density  relates 
the  actual  void  ratio  of  tne  soii  to  the  maximum  and  minimum  void  ratios 
possible  within  the  soil  mass. 

2.1.2  Soil  Benavior  During  Unloading  and  Reloading 

When  an  applied  load  is  removed  from  a  soil  mass,  rebound  will  normally 
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occur  resulting  in  an  increase  in  soil  volume.  The  stress  path  taken  by  the 
soil  mass  during  unloading  will  typically  be  different  from  that  taken  during 
initial  loading.  A  typical  unloading  curve  is  shown  in  Fig.  2.1  as  that 
portion  of  the  curve  lying  between  points  2  and  5.  An  important  fact  in  Fig. 
2.1  is  that  tne  stress  occurring  within  the  soil  mass  is  not  a  single-valued 
function  of  strain.  Instead,  tne  stress  at  a  particular  value  of  strain  may 
be  multi-valued. 

The  term  reloading  refers  to  the  addition  of  a  load  to  the  soil  mass  wmch 
results  in  a  stress  state  previously  experienced  by  the  soil.  A  typical 
reload  curve  is  shown  in  Fig.  2.1  as  that  portion  of  tne  curve  connecting 
points  5  and  3.  When  a  soil  experiences  an  unload-.reload  cycle,  there  will, 
in  general,  be  a  volume  change  associated  with  this  cycle.  As  shown  in  Fig. 
2.1,  tne  unload-reload  cycle  begins  at  point  2  and  ends  at  point  3.  The 
volume  change  which  occurs  during  this  cycle  is  proportional  to  the  difference 
in  the  volumetric  strains  corresponding  to  points  5  and  3  When  the  reload 
path  reaches  point  3  of  Fig.  2.1,  continued  loading  will  follow  a  path  similar 
to  that  for  initial  loading. 

There  will  be  an  energy  loss  associated  with  tne  unload-reload  cycle  as 
apparent  from  the  hystersis  loop  shown  in  Fig.  2.1.  This  demonstrates  the 
effect  of  damping  present  within  a  soil  mass. 

2.1.3  Behavior  of  Soils  in  Simple  Shear 

The  term  simple  shear  means  that  the  soil  is  loaded  in  pure  shear.  The 
behavior  of  soils  when  loaded  in  simple  shear  wili  depend  on  the  initial 
relative  density  of  the  soil.  When  a  soil  of  an  initially  low  relative 
density  is  loaded  in  simple  snear,  a  densif ication  of  the  soil  will  result. 
This  decrease  in  volume  is  due  to  particle  rearrangement,  yielding  and 
fracture.  Densif ication  continues  with  increased  loading  until  a  minimum  void 
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ratio  is  reached.  Continued  loading  will  cause  the  soil  to  fail.  A  dilation 
of  the  soil  mass  is  usually  associated  with  this  failure.  The  dilation  of  tne 
son  mass  occurs  because,  in  order  for  the  soil  to  fail,  grains  must  ride  over 
one  another.  Soils  which  exhibit  such  behavior  are  loose  granular  materials 
and  normally  consolidated  clays.  For  soils  of  initially  high  relative 
density,  little  densif ication  occurs  from  the  application  of  a  simple  snear 
loading.  Rather,  a  dilation  of  the  soil  mass  will  occur  since  tne  void  ratio 
of  the  soil  mass  is  already  near  its  minimum  value.  Soils  whicn  snow  tms 
type  of  behavior  are  dense  sands  and  overconsolidateG  clays.  Typioa*  stress- 
strain  curves  for  different  soils  loaded  in  pure  shear  are  shown  in  Fig.  i.2. 
2.1. A  Soil  Deformation  in  Time 

For  some  soils,  the  total  deformation  resulting  from  the  application  of  a 
load  will  not  occur  instantaneously,  but  over  a  period  of  time.  This  type  of 
deformation  is  referred  to  as  consolidation  and  occurs  in  silts  and  clays. 
Theories  which  predict  the  amount  and  rate  of  consolidation  usually  consider 
tne  soil  to  be  saturated.  The  load  is  initially  transferred  to  the  liquid 
phase  present  in  the  pores  of  the  soil  mass.  This  results  in  an  increase  in 
the  pore  pressure  so  that  steady-state  conditions  in  the  pores  no  longer 
exist.  Over  a  period  of  time,  liquid  will  flow  from  the  pores,  thus  causing  a 
dissipation  of  the  pore  pressures.  This  continues  until  hydrostatic  pressure 
is  acnieved.  As  the  pressure  is  dissipated  from  the  pores,  the  load  will  be 
transferred  to  the  soil  grains  resulting  in  consolidation  of  the  soil  mass. 

The  permeability  of  the  soil  controls  the  rate  at  which  consolidation  takes 
place.  This  behavior  is  termed  primary  consolidation  and  is  shown  in  Fig. 

2.  j. 

Secondary  consolidation  or  creep  is  also  shown  in  Fig.  2.3.  Secondary 
consolidation  is  defined  as  tne  deformation  which  takes  piace  after  the  pore 


Axial  strain  (%) 


Stress-strain  curves  for  loose  and  dense  specimens. 
Medium-fine  sand.  a3  =  30  lb/in.2 :  et  =  0.605  «a  100/o  Dr ; 
e0  =  0.834  ^  20 °o  Dr.  Solid  line,  actual  test  data;  dashed 
li.. .,  extrapolations  based  on  results  of  other  tests.  (After 
Tavlor,  1948.) 


pressures  have  reacned  steady  state  conditions.  Aitnougn  tneones  exist  for 
the  prediction  of  secondary  consolidation,  none  have  yet  found  general 
acceptance. 


2.2  Constitutive  Models  Representing  Soil  Behavior 

Ovci  the  last  two  decades  much  work  nas  been  done  to  develop  constitutive 
models  to  represent  the  load-deformation  behavior  of  soils.  Thus  far,  most  of 
these  models  have  been  used  to  represent  the  behavior  of  dry  or  completely 
saturated  soils.  When  these  models  have  been  used  to  predict  the  behavior  of 
partxy  saturated  soiis,  they  have  yielded  poor  results  for  the  soil  response. 
However,  tnese  models  are  worthy  of  some  attention,  since  they  provide  insight 
to  the  approaches  wmcn  nave  been  taken  to  develop  constitutive  laws 
describing  soil  benavior. 

The  procedure  of  developing  a  constitutive  model  for  soils  has  followed 
two  approacnes.  These  two  approaches  are  termed  phenomenological  and 
micromechamcai  modeling.  The  following  sections  of  this  report  will  discuss 
the  soil  models  obtained  from  these  two  approaches. 

2.2.1  Phenomenological  Models 

Phenomenological  models  are  those  concerned  with  describing  behavior  on 
the  scale  of  tne  experiment.  These  models  treat  soil  as  a  continuum  including 
thousands  to  millions  of  soil  grains  and  pores  in  the  model  representation. 
Phenomenological  methods  or  tneories  include  mathematical  curve  fitting, 
elasticity  theory,  plasticity  tneory  and  visco-elasticity  theory.  The  soil 
constitutive  models  developed  from  these  methods  or  theories  are  discussed 
below  in  detail. 


2. 2. 1.1  Empirical  Models 


A  number  of  models  representing  soil  behavior  have  been  developed  using 
empirical  curve  fitting  methods.  This  approach  entails  making  a  mathematical 
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fit  to  experimental  data.  In  this  manner,  the  response  of  the  soil  due  to 
some  specific  input  may  be  predicted.  Herein,  the  input  and  the  response  are 
those  measured  experimentally. 

Many  researchers  have  taken  the  empirical  approach  to  modeling  the 
pressure-volume  behavior  of  soils.  Herrmann  (1971)  took  such  an  approach  in 
introducing  the  "P-a"  description.  In  his  model,  the  pressure  was  assumed  to 
be  a  function  of  the  specific  volume,  internal  energy  and  the  porosity  of  the 
son.  The  relationship  Herrmann  proposed  is 

v 

p  -  f  (^,  us)  (2.1a) 

v 

a  -  7s  (2.1b) 

m 

wnere  vg  -  specific  volume 

u  »  specific  internal  energy 

v  -  specific  volume  of  the  matrix  material 
m 

In  Eqs.  (2.1),  the  function  f  was  assumed  to  be  identical  to  that  which 
relates  pressure  ana  volume  for  a  mass  composed  entirely  of  matrix  material. 
Carrol  and  Holt  (1972)  proposed  that  it  was  more  reasonable  to  represent  the 
pressure-volume  relationship  for  soils  by 

P  ‘  J  f  (^.  (2’2) 

If  tne  pressure  volume  relationship  for  the  matrix  material  is  known,  the 
problem  reduces  to  determining  the  function 

a  -  q  (P)  (2.3) 

The  determination  of  the  porosity  provides  the  pressure-volume  relationship 
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for  the  complete  material.  This  moaei  used  the  Mie-Gruneisen  equation  of 
state  (see  Carrol  and  Hoit,  1972)  to  describe  tne  function  f.  This  equation 
is 


T 


s 


(2.4) 


where  P  =  pressure 

Pr  »  a  reference  pressure 

u^  =  specific  internal  energy 

u^  =  a  reference  specific  energy 

v  =■  specific  volume 
s 

T  »  Gruneisen  ratio 
g 

The  Gruneisen  ratio  is  given  by 


T«  -  (rc>Cr 


(2.5) 


where  T  -  temperature 

-  specific  heat  at  constant  volume 

A  polynominal  fit  was  then  used  to  describe  the  function  g.  For  situations 
where  the  variation  m  internal  energy  is  less  important,  Butkovicn  (1973) 
developed  a  model  relating  the  porosity  to  the  applied  pressure.  The 
expression  ne  obtained  is 


l 

a. 


,  ,,  1  P 

'  ’  u  ’  5")  Vr> 

_ o _ c__ 

Vr' 


(2.6) 


where  a  -  initial  value  of 
o 

Pg  -  the  pressure  required  for  the  onset  of  pore  closure 


*  the  pressure  requirea  for  pressures  within  the  range 
The  Eq.  (£.4)  is  applicable  for  pressures  witnin  the  range  between  Pg  and  P  . 
For  pressures  less  tnan  P  ,  the  pressure-volume  relationship  for  the  matrix 
material  is  used.  In  Butkovicn’s  work,  the  pressure-volume  relationship  for 
the  matrix  material  was  assumed  to  be  given  by  soil  unloading  data.  A 
polynomial  fit  to  initial  loading  data  was  used  to  determine  the  pressure 
volume  relationship  for  pressures  lying  between  Pg  and  P  .  Otner  empirical 
models  describing  the  pressure-volume  relationship  of  soils  have  been 
developed,  but  the  models  cited  above  are  representative  of  this  work. 

Other  models  have  been  developed  wmch  made  mathematical  fits  to  deviator 
stress-strain  data.  The  simplest  model  of  tms  type  is  obtained  by 
approximating  deviator  stress  versus  strain  data  by  a  series  of  piecewise 
linear  curves.  This  type  of  approximation  is  shown  in  Fig.  2.4.  More 
sophisticated  mathematical  fits  such  as  hyperbolas  and  cubic  splices  have  been 
used  to  relate  tne  deviator  stress  to  the  strain  data.  The  most  popular  of 
these  methods  is  the  hyperbolic  stress-strain  model  developed  by  Duncan  and 
Chang  (1970).  This  model  is  based  on  the  finding  that  the  deviator  stress 
versus  strain  curves  for  a  number  of  soils  can  be  approximated  by  hyperbolas 
like  the  one  shown  in  Fig.  2.5.  This  hyperbola  may  be  represented  by 
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~  °3  ^ULT 


(2.7) 


wnere  n  »  deviator  stress 

e  *  strain 

E  -  initial  tangent  modulus 

The  Eq.  (2.7)  may  be  transformed  so  that  it  will  plot  as  a  straight  line. 


This  transformed  equation  is 
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A  plot  of  Eq.  (2.8)  is  shown  in  rig.  2.6.  Other  empirical  relations  are  used 
to  account  for  the  variation  in  soil  strength  with  confining  pressure  and 
modulus  values  for  loading  and  unloading.  Because  of  these  relationships,  the 
hyperbolic  model  requires  a  number  of  parameters  for  its  use.  The  hypobolic 
model  also  fails  to  realistically  model  actual  soil  behavior  at  or  near 
failure. 

In  general,  there  are  some  basic  problems  associated  with  soil  models 
developed  from  empirical  models.  First,  an  empirical  model  cannot  be  expected 
to  provide  reasonable  predictions  of  soil  behavior  when  the  soil  and  site 
conditions  being  modeled  deviate  greatly  from  those  used  to  calibrate  the 
model.  Second,  this  type  of  model  cannot  oe  expected  to  provide  any  insight 
to  the  actual  pnysical  deformation  mechanisms  acting  within  the  soil  mass. 
Despite  these  shortcomings  empirical  models  are  frequently  used  due  to  their 
simplicity. 

2. 2. 1.2  Nonlinear  Elastic  Models 

Some  models  representing  soil  behavior  have  been  developed  using  nonlinear 
elasticity  theories.  However,  these  theories  have  not  found  widespread  use 
since  tneir  predictions  of  unload  behavior  do  not  represent  actual  soil 
behavior.  For  cases  where  initial  loading  is  of  interest,  nonlinear 
elasticity  theories  may  provide  reasonable  predictions  of  soil  response. 

Hyperelastic  constitutive  laws  have  been  used  to  represent  soil  behavior. 
These  models  use  constitutive  laws  obtained  by  the  differentiation  of  a  strain 
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energy  function.  Different  orders  of  hyperelastic  models  are  obtained  by 
retaining  the  higher  order  derivatives  obtained  from  the  strain  energy 
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Truesdale  nas  proposed  a  rate  theory  which  states  that  the  rate  of  change 
of  stress  is  a  function  of  the  rate  of  change  of  strain.  This  is  known  as  tne 
hypoeiastic  formulation.  At  present,  this  formulation  has  not  been  used  in 
representing  the  load-deformation  behavior  of  soils. 

2. 2. 1.3  Elastic-Plastic  Models 

Elastic-plastic  theory  had  been  widely  used  in  soil  modeling.  Recently, 
many  constitutive  models  for  soils  have  been  presented  which  use  this  theory. 

In  general,  this  type  of  model  assumed  a  yield  criterion  of  the  form 
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F  (Ty.  .y.  X) 


(2.9) 
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where  T^.  -  stress  tensor 

-  plastic  strain  tensor 
X  -  work  hardening  parameter 

p 

When  the  above  equation  is  not  satisfied  (F  »  X)  <  0)  the  material 

is  said  to  behave  elastically.  When  Eq.  (2.9)  is  satisfied,  the  behavior  is 
said  to  be  elastic-plastic.  Furtner  deformation  beyond  the  point  where  Eq. 
(2.9)  is  satisfied  occurs  at  a  combination  of  elastic  and  plastic  strains, 
prescribed  by  an  assumed  flow  rule.  The  yield  surface  is  typically  described 
in  stress  space  as  shown  in  Fig.  2.7.  For  a  known  stress  point  inside  this 
region,  tne  strains  are  found  using  elastic  constitutive  laws.  When  the 
stress  point  lies  on  the  yield  surface,  the  total  strain  is  a  combination  of 
elastic  and  plastic  strains.  For  a  stress  point  lying  on  the  yield  surface, 
further  loading  may  cause  the  surface  to  expand,  translate  or  both,  according 
to  the  work-hardening  rule  assumed.  Unloading  may  be  elastic  or  elastic-; 
plastic.  An  example  of  the  movement  of  the  yield  surface  in  the  principal 
stress  space  is  shown  in  Fig.  2.8.  The  material  is  initially  unstressed  at 
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point  0  and  is  then  loaded  to  point  1 .  The  benavior  of  the  material  when 

following  this  load  path  is  elastic,  since  it  lies  within  tne  initial  yieid 

surface.  At  point  1 ,  the  stress  point  lies  on  the  yield  surface  such  that 

further  loading  will  result  in  behavior  that  may  be  elastic-plastic  or 

elastic,  depending  on  the  stress  path  taken.  As  the  material  is  loaded  along 

the  stress  path  connecting  points  1  and  2  in  Fig.  2.8,  the  material  will 

exhibit  elasticrpiastic  behavior.  The  stress  point  along  this  path  remains  on 

the  yield  surface,  with  the  yield  surface  expanding,  translating  or  botn.  In 

Fig.  2.8,  the  tensor  a  .  would  be  non-zero  and  no  expansion  of  the  yield 

i  J 

surface  would  occur.  Combined  hardening,  in  which  the  yield  surface  may 
translate  as  well  as  expand,  is  snown  in  Fig.  2.8.  Further  loading  of  the 
material  from  point  2  to  point  3  of  Fig.  2.8  will  result  in  elastic  behavior 
because  the  load  path  taken  lies  within  the  new  yield  surface.  Continued 
loading  from  point  3  will  result  in  elastic  or  elastic-plastic  behavior, 
depending  on  the  load  path  taken. 

Schofieid  ana  Wroth  (1968)  developed  an  elastic-plastic  soil  model,  known 
as  the  "cam-clay"  model  which  accounts  for  the  volume  deformation  and  strain¬ 
hardening  of  soils.  Tne  basis  of  their  model  is  an  incremental  flow  rule 
whicn  balances  the  irreversible  work  occurring  during  deformation  against  a 
mecnanism  for  the  frictional  loss.  Their  flow  rule  is  given  by 

P  aV  P  P 

Q  -  (Atg)  -  P  (— )  ,  UP  (Ats)  (2.10) 

wnere  V  =*  volume 

P  -  pressure 

Q  -  measure  of  shear  stress 
t  -  measure  of  shear  strain 
u  -  friction  parameter 


Here  tne  superscript  denotes  plastic  or  elastic  portions  of  tne  quantities 
indicated.  The  elastic  volume  deformation  during  hydrostatic  deformation  is 
given  by 


,V  AP  V  6  A P 

d(— )  =  -  A  —  A  (— )  =  -  A  —  (2.11) 

m  m 


wnere  V  =  volume  of  matrix  material 

m 

A  =  a  constant 

As  yielding  is  occurring,  the  total  volume  change  is  given  by 


„  ,V  AP  V  AP 

u  y  )  =  ~  B  p  A  (y)  =  B  p  (2.12) 

mm 

where  B  =  a  constant 

Tne  assumption  of  an  associated  flow  rule  gives  the  following  equation  for  the 
yield  surface: 


In  Eq.  (2.13)  P*  is  the  intercept  of  the  yield  surface  with  the  P  axis  as 
shown  in  Fig.  2.9.  An  important  assumption  of  the  "caimclay"  model  is  that 
tne  plastic  volume  deformation  for  non-hydrostatic  stress  states  is  the  same 
as  for  hydrostatic  states  but  with  the  p  replaced  by  P*.  Thus,  tne  plastic 
volume  deformation  is  given  by 


,  V  AP* 

A(— )  =  -  (B  -  A)  “pqp  (2.14) 

m 

Tne  Eqs.  (2.10),  (2.11),  and  (2.12)  form  a  system  of  equations  from  which 
strain  increments  may  be  determined  from  stress  increments,  or  vice-versa. 
The  constants  B  and  A  are  determined  experimentally.  Hydrostatic  loading 
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corresponas  to  a  movement  along  the  P  axis  shown  in  Fig.  2.9.  */hen  yielding 
occurs,  a  non-nyaroelastic  loading  will  cause  tne  yield  surface  to  change  in 
accordance  with  Eq.  (.2.10).  Movement  of  tne  yield  surface  is  shown  in  Fig. 
2.9.  The  critical  state  line  is  snown  m  Fig.  2.9  as  tne  line  connecting 
points  of  zero  slope  for  all  possible  yield  surfaces.  This  separates  yielding 
into  densif ication  and  dilation.  Densif ication  with  strain-hardening  occurs 
to  tne  right  of  the  critical  state  line,  while  dilation  witn  strain  softening 
occurs  to  the  left. 

Tne  "cam-cxay"  model  nas  proved  useful  in  representing  soil  behavior.  In 
this  model,  however,  elastic  shear  stresses  and  soil  cohesion  are  completely 
neglected.  '  The  assumption  of  an  associated  flow  rule  is  also  made.  This 
assumption  gives  a  plastic  strain  vector  normal  to  the  yield  surface. 
Contributions  sucn  as  those  by  Mandl  and  Luque  (1970)  and  Frydman  et  al. 

(1973)  nave  shown  that  normality  of  plastic  flow  is  neither  a  mathematical 
necessity  nor  supported  by  experimental  evidence.  The  "cam-clay"  model 
predicts  no  non-recoverable  deformations  under  hydrostatic  loadings.  This  is 
not  representati ve  of  soils.  Unloading  is  elastic,  wmcn  is  not  descriptive 
of  actual  soil  benavior. 

Sandier  and  Baron  have  introduced  the  "cap"  model  to  describe  the  benavior 
of  sons.  This  model  is  based  on  the  classical  plasticity  model,  defined  by  a 
yield  surface  and  a  strain  rate  vector.  Inspection  of  this  yield  surface  shows 
that  tnree  modes  of  soil  behavior  are  possible.  These  being  elastic,  failure, 
and  cap  behavior.  Elastic  behavior  occurs  when  the  stress  point  lies  in  the 
region  contained  by  tne  stress  coordinate  axes,  the  failure  envelope,  and  the 
cap  surface.  The  benavior  in  this  region  is  considered  to  be  linearly 
elastic.  The  failure  model  of  behavior  occurs  when  tne  stress  point  lies  on 
the  failure  envelope.  This  failure  envelope  is  assumed  to  be  fi-xed  and  is 
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given  by 


4  =  A  -  Ce'3BP 


(2,14) 


where  =  second  invariant  of  the  deviatonc  stress  tensor 


A , B, C  =  material  constants 

The  moael  assumes  an  associated  flow  rule,  thus,  the  plastic  strain  during  the 
failure  mode  of  behavior  is  composed  of  a  shear  component  and  a  dilatant 
component.  The  cap  mode  of  behavior  occurs  when  tne  stress  point  lies  on  the 
cap  surface  and  continued  loading  results  in  an  outward  movement  of  the  cap. 
The  motion  of  tne  cap  is  related  to  the  plastic  strain  by  a  hardening  rule. 

The  equation  for  the  cap  surface  is 


(P  "  Pa)  +  \  D^I ^ J 2 


<Pb  "  Pa> 


(2.15) 


where  P  ,P  =*  pressures  corresponding  to  points  a  and  b  as 

ct  D 


on  the  Yield  surface 

P  =  mean  hydrostatic  stress 


1,  -  trace  of  the  stress  tensor 


=  second  invarient  of  the  deviatonc  stress  tensor 


D  =  a  constant 


Tne  position  of  tne  cap  is  defined  by  specifying  one  of  the  quantities,  P  or 

cl 


P^.  Tne  cap  is  related  to  the  strain  history  of  tne  material  through  a  strain 


hardening  rule  given  by 


-P 

e 

y 


-3NP- 

N  (l  -  e  B) 


(2.16) 


where  M,N  =  material  constants 
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Wnen  tne  stress  pc;. it  lies  on  eitner  tne  failure  envelope  or  tne  cap  surface, 

-P 

tne  value  of  e  cnanges  exactly  as  tne  plastic  volumetric  strain  for  a 

v 

stress  point  on  the  cap  surface,  tne  plastic  strain  rate  vector  will  oe 

directed  as  snown  in  Fig.  (2.o).  The  position  of  the  plastic  strain  rate 

vector  implies  tnat  it  consists  of  an  irreversible  decrease  in  volume  in 

conjunction  witn  an  irreversible  snear  strain.  This  decrease  in  volume 

represents  volumetric  nysteresis  observed  in  soil  during  compaction,  As  the 
cap  moves  forward,  tne  compaction  resulting  from  tne 

„ ,  -p 

associated  fiow  will  lead  to  an  increase  in  tne  Cap  parameter  .  By 

v 

Eq.  (2.1o)  this  leads  to  an  increase  in  P  ,  resulting  in  a  movement  of  tne  cap 
to  tne  rignt .  when  the  stress  point  lies  on  the  failure  surface,  tne  plastic 
strain  rate  vector  will  be  directed  upwards  and  to  tne  left  as  shown  in  Fig. 
2.6.  The  plastic  strain  rate  vector  indicates  an  increase  in  volume 
associated  witn  tne  movement  along  tne  failure  surface.  The  dilatancy  will 
leaa  to  a  decrease  in  tne  cap 

-P 

parameter,  resuxting  in  a  leftward  movement  of  tne  cap  by  Eq. 

y 

(2.1o).  Tne  backward  movement  of  the  cap  is  limited  by  tne  point  where  it 
intersects  tne  stress  point  lying  on  tne  failure  surface. 

Tne  basic  cap  model  described  has  been  modified  to  include  viscous  dampin„ 
ana  strain  naraening.  Tne  viscous  cap  model  is  used  to  represent  materials 
wmch  exmbit  hysteresis  during  cyclic  loading.  This  model  was  formulated  by 
introducing  linear  viscous  damping  into  tne  elastic  portion  of  tne  cap  model. 
The  parameters  wmch  define  the  non-plastic  portion  of  tne  model  are  an 
instantaneous  modulus  G^,  a  long  term  modulus  Gg,  and  a  relaxation  rate  x.  The 
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tnree  parameters  are  related  through  the  equation 


3  gf+gv 


(2.17) 


(Gp 
a  F 


(2.18) 


wnere  Gy,  Gp-  spring  moduli 

-  long  term  modulus 
G^  -  damping  constant 
t  -  relaxation  rate 

The  deviatoric  stress-strain  relation  for  tne  viscous  cap  model  is  given  by 


d(d  V  )  2  G  dV .  -  S 

U  s  11  1. 


JF  dt 


(2.19) 


wnere  S  -  deviatoric  stress  tensor 
ij 

dv  -  viscoelastic  deviatoric  strain  tensor 
ij 

To  determine  the  parameters  Gp,  Gg,  and  t, cyclic  triaxiai  data  are  used.  A 
Kinematically  hardening  failure  envelope  has  been  added  to  the  general  cap 
model  by  replacing  the  stress  tensor  x^  by  ( x^  -  a-  ) .  Here  is  a  tensor 
wnose  components  are  memory  parameters  defining  the  translation  of  the  failure 
surface  in  stress  space.  In  the  model,  it  is  assumed  that  kinematic  hardening 
occurs  only  in  shear,  yielding  tne  relation 


G  -  0 

kK 


(2.20) 


'he  kinematic  naroenmg  rule  which  governs  tne  memory  parameters  G^  is  of  the 


form 
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wnere 


a  -  -  tensor  defining  tne  translation  of  tne  failure  surface 
ij 


ij 


ij 


stress  tensor 

a  work  naroening  parameter 
strain  tensor 


d,  -  aiviatoric  plastic  strain  tensor 

KX 

In  oraer  to  represent  tne  nonlinear  penavior  of  sens  at  or  near  failure,  it 
is  necessary  to  assume  a  nonlinear  hardening  rule.  A  simple  rule  of  this  type 
wmen  gives  reasonable  behavior  at  all  stress  levels  is  given  as 


U 


C  F  e  . 
a  a  lj 


(2.22) 


wnere 


C  -  a  constant 

a 


1  -  (  t  .  ~  a  )  a 

F  -  maximum  [0,- - ^ ^ - -*— ] 

a  2N  (/J,_  -  N ) 


2F 


/J-F-  A  -  Ce 


-3BP 


Ny  -  a  constant  defining  tne  size  of  tne  yield  surface 
Here  F„  is  related  to  tne  proximity  of  tne  yield  surface  to  the  failure 

U 

surface,  and  tne  location  of  tne  stress  point  on  the  yield  surface.  For  a 


ij 


j,  F  wnl  be  equal  to  1.0.  Therefore,  from  Eq.  (2.22),  it  is  found  tnat  Ca  is 
the  inelastic  slope  for  the  initial  yielding  of  tne  material  in  shear. 
will  decrease  for  continued  yield,  and  is  equal  to  zero  when  the  stress  point 
reaches  tne  failure  surface.  Upon  unloading  from  tne  failure  surface,  tne 


vV-' 


va iue  of  F  win  increase,  reacrung  a  value  of  2  upon  reyieiding.  Finally,  one 
cap  raoaex  nas  been  modified  to  represent  tne  behavior  of  saturated  soils  using 
tne  effective  stress  approach.  This  modification  is  straightforward  ana  is 
achieved  by  replacing  tne  stress  tensor  by 


T*  -  T  -  Ub 

IJ  1J  U 


(2.23) 


where 


effective  stress  tensor 


x  ,  -  total  stress  tensor 
ij 

u  -  pore  pressure 


d  Kronecker  -  [’  J } 
ij  l0,  1  * 


The  cap  monel  has  been  used  successfully  to  moael  several  soils.  However, 
there  are  some  difficulties  associated  with  it.  A  major  problem  is  that  a 
large  number  of  parameters  must  be  determined  from  experimental  data  and  their 
determination  may  require  special  tests.  Another  problem  is  the  assumption  of 
an  associated  flow  rule.  This  assumption  is  not  necess'anly  correct  for 
soils. 

Other  elastic-plastic  constitutive  models  for  soils  have  been  developed. 
These  models  use  different  yield  surfaces  as  a  nonnassociated  flow  rule. 
However,  the  methodology  used  to  formulate  these  models  is  the  same  as  for 
those  already  described.  The  problem  with  classical  plasticity  theory  is  that 
tne  predicted  response  of  a  system  is  rate  independent.  It  has  been 
established  that  the  response  of  a  soil  is  rate  dependent. 
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2. 2. 1.4  Viscoelastic  Models 

Viscoelastic  moaej-s  aescrioing  soix  oenavior  nave  not  appeared  in  tne 
-literature  as  mucn  as  those  formulated  using  elastic-plastic  tneory.  However, 
eiastic-piastic  tneory  is,  in  fact,  a  special  case  of  viscoelastic  tneory. 
Using  tne  more  general  material  model  provided  from  viscoelastic  tneory,  soil 
damping  ano  rate  dependence  may  be  accounted  for.  As  seen  for  the  "cap 
mooei",  results  improved  wnen  viscous  damping  was  introduced  into  tne  model. 
2.2.2  Micromecnanical  Models 

Mecnamstic  modeling  of  soils  nas  been  approacned  from  two  different 
viewpoints.  One  apphoacn  nas  been  to  treat  tne  soil  as  an  assemoiage  of 
particles  in  contact.  The  particles  within  a  soil  mass  may  be  random  in  shape 
and  size;  tnerefore,  to  use  this  approacn  some  assumptions  as  to  size  and 
snape  must  usually  be  made.  Once  a  model  representing  the  soil  mass  has  been 
cnosen,  the  solution  consists  of  representing  the  deformed  geometry  of  the 
particles  in  contact.  The  other  approach  to  mechanistic  modeling  has  been  to 
consider  the  3oii  as  composed  of  a  matrix  material  containing  voids.  A 
solution  to  this  problem  consists  of  modeling  the  deformation  of  tne  voids 
contained  in  the  matrix  material. 

Mecnamstic  models  have  been  formulated  on  two  scales.  One  scale  has  been 
intermediate  to  that  of  tne  experiment  and  the  grains  and  the  pores  within  the 
soil  mass,  while  this  scale  may  be  very  small  compared  to  the  scaie  of  tne 
experiment,  it  may  be  quite  large  in  comparison  to  the  size  scale  of  the 
grains  ano  pores.  On  tms  scale,  the  benavior  observed  may  be  that  of  many 
grains  and  pores  and  may  best  be  described  by  the  use  of  a  phenomonological 
tneory.  The  other  scaie  wmch  is  used  in  mechanistic  modeling  is  termed  tne 
micro-scale.  On  this  level,  models  are  formulated  at  the  scale  of  the  grains 
and  pores  and  are  concerned  witn  describing  the  actual  deformation  mecnamsms 
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present,  on  crus  scaxe. 

Tne  void  deformation  models  nave  oeen  formulated  on  ootn  tne  intermediate 
scaxe  and  micronscaxe.  Modeling  of  objects  in  contact  nas  usually  been  done 
on  tne  micro-scale. 

2.2.2. 1  Contact  Mooeis 

«’hen  a  mass  composed  of  a  number  of  particles  in  contact  is  subject  to  an 
externally  applied  load,  tne  deformation  resulting  from  tne  load  is  due  to 
grain  movement  and  grain  deformation.  The  movement  of  tne  grains  will  be 
controlled  by  mterparticie  friction,  cohesion  between  adjacent  particles,  and 
tne  initial  porosity  of  tne  mass.  The  grain  deformation  will  be  greatest  at 
areas  of  contact  between  adjacent  grains,  and  may  be  elastic  or  elastic-! 
pxastic,  depending  on  the  stress  level  present  in  the  grains.  In  addition, 
tne  grains  may  fracture,  thus  changing  the  number  and  tne  snape  of  the  grains 
and  increasing  the  number  of  contacts. 

Models  used  to  describe  this  behavior  usually  consider  the  soil  grains  to  be 
spnencal  in  snape.  The  load-deformation  behavior  of  the  spneres  themselves 
is  considered  to  be  that  of  an  elastic  material.  Further  simplifications  are 
obtained  by  neglecting  friction,  cohesion,  and  tangential  forces  acting  on  tne 
contacts  between  grains.  With  these  simplifications,  a  logical  step  is  to  use 
Hertzian  contact  theory,  by  which  the  movement  of  adjacent  spheres  relative  to 
one  anotner  may  be  determined.  Two  spheres  in  contact  are  shown  in  Fig.  2.10. 
The  z  axis  is  positioned  at  tne  centerline  of  the  contact.  The  soiid  lines 
represent  tne  deformed  configuration  of  tne  spneres,  wmie  tne  dashed  xines 
represent  the  undeformed  spneres.  From  Hertzian  contact  theory,  tne 
deformation  axong  tne  centerline  of  contact  for  each  sphere  is  given  by 
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(2.24b) 


wnere 


1  2 

u  ,  u  -  aeformation  along  the  contact  centerline  for 

2  Z 


spneres  1  ana  2  respectively 


vl*v2  *  Poisson's  ratio  for  spheres  1  ana  2  respectively 


-  eiastic  modulus  for  spneres  1  ana  2  respectively 


a  »  raaius  of  tne  contact  surface 


F  -  force  transmitted  across  the  contact 


The  aeformations  u1  and  ^  are  shown  m  Fig.  2.10.  The  raaius  of  tne  contact 


area  between  the  spheres  is  given  by 


3tt  P  RjR2  (1  ~  vj2) 
4  ^  +  R2  E1 


(2.25) 


where  -  the  radii  of  spneres  1  and  2  respectively. 


Using  Eqs.  (2.24)  and  (2.25),  the  deformation  of  an  assemblage  of  spheres  may 


be  determinea  when  the  force  transmitted  across  each  contact  is  known.  Ko  and 


Scott  (1967)  have  solved  this  problem  for  the  case  of  an  assmebly  of  spheres 


in  ideal  packing  configurations,  under  conaitions  of  hyarostatic  loading.  All 


tne  spneres  were  considered  to  be  of  equal  radii  and  of  the  same  material 


properties.  The  solution  is  given  by 


A V  ,  3  3C  (1  -  vZ)  P 
A  -  E 

•  « 


(2.26) 


wnere  V  -  volume  of  soil  mass 


C  -  1  ,  for  sc  (3imple  cubic)  packing 


/  2/4,  for  fee  (face  centered  cubic) 

P  -  nyarostatic  pressure 
E  -  elastic  modulus 
v  -  Poisson's  ratio 

As  seen  from  Eq.  (2.26),  the  term  C  accounts  for  tne  initial  density  of  tne 
mass,  giving  smaller  volumetric  strains  for  the  denser  packing  configurations. 

However,  this  model  predicts  larger  tnan  actual  deformations  for  the  simple 
cuDic  (sc)  conf iguration,  while  predicting  smaller  than  actual  deformations 
for  tne  face  centered  cubic  (fee)  packing  configuration.  To  correct  this,  Ko 
ana  Scott  used  a  combination  of  sc  and  fee  blocks  to  represent  the  initial 
porosity  of  tne  soil.  By  assuming  a  distribution  of  gram  contact  pressures 
and  an  effective  contact  radius,  pressure-volume  relationships  for  sands  of 
three  initial  porosities  were  generated.  The  results  obtained  are  shown  in 
Figs.  2.11  and  2.12  along  with  tne  limiting  cases  of  sc  and  fee  packing 
configurations.  A  major  shortcoming  of  Ko  and  Scott's  model  is  that  the  path 
the  soil  takes  during  unloading  is  the  same  as  that  for  loading,  which  is  not 
representati ve  of  actual  soil  behavior.  Warren  and  Anderson  (1973)  have 
formulated  a  contact  model  in  which  initially  some  of  the  spheres  are  not  in 
contact.  The  pressure-volume  relationship  obtained  is  given  by 


AV 

V 


3 


N 

c 


li_L 


1  2/3 


(2.27) 


wnere  -  number  of  grains  in  a  typical  cross-section 

M  ■  The  number  of  contacts  transmitting  force  across  the 
c 

typical  cross-section 

As  loading  progresses,  more  grains  come  into  contact,  until  at  some  critical 


pressure,  all  grains  make  contact.  It  is  apparent  from  Eq.  (2.27)  that  as  the 
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number  of  contacts  is  increased,  the  amount  of  volume  deformation  resulting 


from  an  increase  in  pressure  win  decrease.  The  mooel  will  predict  unloading 


a^ong  a  patn  different  from  tnat  of  loaning,  as  long  as  the  grains  were  not 


all  initially  in  contact.  Tne  difficulty  witn  this  model  is  the  determination 


of  tne  vaiue  of  N  .  The  variation  m  the  value  of  N  wnicn  occurs  during 

C*  '•rf 


loading  corresponds  to  the  rigid  body  motion  of  tne  particles  within  the  soil 


mass.  This  model  does  not  attempt  to  describe  the  actual  grain  motion  within 


the  soil  mass,  but  ratner  tne  parameter  N  is  cnosen  to  fit  experimental  data. 


Some  mooeis  of  granular  media  include  friction  of  the  contact  between 


grains..  Rowe  (1962)  nas  considered  tne  shearing  of  various  assemblages  of 


spheres.  Using  a  minimum  energy  criterion,  he  arrived  at  the  stressndilatancy 


equation 


31  2  *u  1  +  V  ‘ 

—  -  tan  (45°  -  —  ) 


(2.28) 


wnere  a.  =*  maximum  principal  stress 


z  »  j.east  principal  stress 


■  undramed  angxe  of  shearing  resistance 


ue  -  the  axial  strain  increment 


This  equation  holds  only  for  the  case  wnen  the  intermediate  principal  stress 


is  equal  to  the  least  principal  stress.  Rowe  states  tnat  tne  angle  $  ,  must 


be  replaced  by  an  effective  angle  of  shearing  resistance  q  ,.  in  order  to  matcn 


experimental  data.  Test  conditions  may  be  created  so  that  many  values  of  tne 


undramed  angle  of  shearing  resistance,  <p  ,  may  be  obtained  for  the  same  soil 


sample.  However,  with  pure  pressure  measurements  during  the  test,  tne  value 


of  the  effective  angle  of  shearing  resistance  may  be  determined  and  this  value 


nas  been  found  not  to  vary  witn  test  conditions.  The  Eq.  (a. 28)  does  not 
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account  for  compaction  during  non-nyarostatic  loading.  Barden  et.  ai.  (1969) 
used  Eq.  (2.28)  to  formulate  a  plastic  flow  rule  and  a  set  of  yield  surfaces. 
Tney  tested  tne  penavior  of  sane  m  plane  strain  ano  found  tnat  tne  yield 
criterion  and  plastic  potential  did  not  coincide.  This  implies  non-normauity 
of  fxow.  However,  it  was  found  that  tne  volumetric  strain  was  suitably 
predicted  oy  tms  model.  Nemat-Naser  formulated  a  model  to  represent  the 
benavior  of  granuxar  material  undergoing  snear  ioaemg.  This  model  represents 
dilation  and  densif ication  which  occurs  curing  shear.  This  is  done  by 
defining  the  dilatancy_angle  6,  wruen  defines  tne  position  of  a  microscopic 
shear  plane  with  respect  to  the  observable  macroscopic  snear  plane.  In  this 
model,  it  is  assumed  that  the  actual  shearing  takes  piace  on  many  microscopic 
snear  planes  rather  than  on  one  macroscopic  shear  plane.  To  formulate  the 
modei,  Nemat-Naser  considered  a  sample  of  soil  for  whicn  failure  takes  place 
along  one  microscopic  shear  plane  S'-S'.  Summing  forces  on  the  plane  S'-S' 
gives  tne  following  equations 


F* 


T  tan  <p.^  cos  (<p  -  v) 
sin  ip 


t.  dfl  (p 


a 


(8  -  v) 


(2.29) 


where  T  »  total  snear  force  on  macroscopic  sample 
6  -  dilatancy  angle 

F*  »  frictional  force  along  microscopic  shear  plane 
Tne  angles  <p  and  ■?  are  defined  by  tne  equations 

t  -  a  tan  6  (2.30a) 

x*  -  a*  tan  (2.30b) 


where 


t 


snear  stress  acting  on  macroscopic  snear  plane 


s 


a  =  normal  stress  acting  on  macroscopic  snear  plane 
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x*  =  snear  stress  acting  on  microscopic  snear  plane 


o*  =  normal  stress  acting  on  microscopic  snear  plane 


The  Eq.  (2.29)  was  derived  by  considering  all  stresses  in  Eqs.  (2.30)  to  be 


acting  on  tne  same  element  of  area  along  tne  microscopic  snear  plane.  3y 


considering  tne  rate  of  energy  dissipation  whicn  occurs  as  the  block  slides 


along  the  plane  S'-S',  Neraat-Naser  obtained  tne  equations 


ian  ?  cos  (f  -  v)  ^ 


(2.31  ) 


sm  q  sin  v 


w  »  rate  of  energy  dissipation 
wnere  V  »  tne  volume  of  tne  macroscopic  sample. 


Tne  dot  denotes  tne  time  derivatives.  The  following  approximations  are  made 


concerning  w: 


w  -  w'  +  w" 


(2.32a) 


•w'  -  xY 


(2.32b) 


-!  -  V 

sm  $  cos 


(2.320 


wnere  Y  «  the  rate  of  snear  deformation  on  the  macroscopic  sample. 


Use  of  these  approximations  in  Eq.  (3.31)  gives  the  equation 


i  vi 


cos  (<J>  +  v  )  sin  v . 

_  u  _ i _ i 


Vi  * 


cos  <P 


(2.33) 


Eq.  (2.33)  applies  to  one  microscopic  failure  pxane.  If  one  microscopic 


snear  plane  is  denoted  by  l,  Eq.  (2.33)  is  written  as 


1  Vi 


cos  ($  +  v . )  sin  v . 
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tne  volume  fraction,  V  ,  of  tne  family  of  particles  naving  a  ailatancy  angle, 
is  defined  by, 


p±(v)  -  v 


(2.35) 


wnere  P  -  volume  fraction  of  family  of  particles  having  ailatancy 
angle 

Tne  restriction  if  P  is  the  following 


J  Pi 


(v)  dv  ■  1 


(3.36) 


In  Eq.  (3-36)  v  ♦  ,nd  v.-  form  the  range  of  variation  of  the  dilatancy  angle, 

o  dUU  1 

V.  Using  Eqs.  (2.35)  and  2.36)  in  Eq.  (2.34),  Nemat^Naser  arrived  at  the 


final  result  given  by 


11m  _ j 

V  y  cos 


o 

I  tp  J'  coa  (<Pu  *  v)  sin  V  dv. 


(3-37) 


The  Eq.  3-37  contains  ail  experimentally  observed  behavior  of  granular 
material  in  simple  shear.  However,  the  accuracy  of  the  predictions  made  by 
Eq.  (3.37)  will  depend  on  the  chosen  form  of  the  distribution  function  P(V). 
This  distribution  function  may  be  very  difficult  to  determine  for  an  actual 
soil  sample.  Anotner  shortcoming  cp  this  model  is  tnat  the  individual 
particles  within  the  sample  mass  are  considered  to  be  rigid.  Wilkins  (1970) 
tooK  a  different  approach  to  develop  a  theory  for  tne  snear  strength  of  a 
granular  media.  He  used  an  empirical  curve  fitting  method  and  Rowe's  Eq. 


) 

; ^ . ^3 ;  co  predict  tne  number  of  unstable  contacts  in  a  granuiar  assemblage  as  < 
a  function  of  tne  stress  ratio.  According  to  tms  approacn,  wnen  all  tne 

I 

contacts  on  a  grain  become  unstable,  tne  grain  is  no  xonger  considered  to 
contribute  to  the  system  and  it  effectively  becomes  a  void,  When  tne  number 
of  voids  not  supporting  any  stresses  is  equal  to  the  number  of  particles  which 
continue  to  carry  loads,  the  medium  is  assumed  to  fail.  Altnough  this  attempt 
is  interesting,  it  becomes  unattractive  due  to  its  empirical  nature.  Volume 
cnanges  and  stress-strain  relations  are  neglected  in  Wilkins'  formulation. 

Other  contact  models  have  been  developed  for  which  the  plastic  flow  of  the 
bodies  in  contact  are  considered  to  be  important.  Kaxar  and  Chaklader  (1967) 
nave  solved  this  problem  for  spheres  in  a  variety  of  packing  configurations. 

In  this  model,  it  is  assumed  tnat  the  particle  surfaces  whicn  are  not  in 
contact  remain  spherical.  They  solved  this  for  a  simple  cubic  packing. 

Assumptions  were  that  tne  volume  of  the  spheres  remain  constant,  that  the 
contacts  transmit  the  load  applied  to  the  assembly,  and  that  the  material  near 
tne  contact  is  in  a  state  of  uniaxial  stress.  The  material  of  the  contacts 
was  allowed  to  yield  until  the  stress  developed  at  the  contacts  was  balanced 
by  the  applied  pressure.  The  relationship  tnat  Kakar  ano  Chaklaoer  obtained 
is  given  by 
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(2.38) 


wnere  ?  -  applied  pressure 

Y  »  yield  strengtn 

V  »  volume  of  tne  assembly 


Tne  £q.  (3-38)  is  valid  until  tne  contact  areas  toucn,  tnus  forming  a  new 
geometry.  The  results  snown  are  for  a  simple  cubic  packing  configuration  ano 
tnese  snow  tnat  tne  yielding  model  predicts  larger  strains  for  a  given  loan 
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tnan  tnose  contained  from  the  elastic  Hertzian  contact  model.  The  actual 
stress-strain  curve  will  iikeiy  fall  in  between  tnose  by  the  Hertzian  and 
Kakar  and  Chaklader  theories.  Not  all  points  within  a  sphere  will  yield  at 
once;  tnus ,  the  actual  behavior  is  stiff er  than  that  predicted  by  the  complete 
yielding  model  as  formulated  by  Kakar  and  Chaklader. 

2. 2. 2. 2  Void  Deformation  Models 

One  approacn  to  modeling  soils  has  been  to  consider  the  soil  as  a  mass 
composed  of  a  matrix  material  and  voids.  The  deformation  resulting  from  the 
application  of  loads  to  a  material  of  this  type  will  aepend  on  the  materials 
making  up  the  matrix  and  voids,  the  size  and  the  shape  of  the  voids,  and  the 
column  fraction  of  the  voids.  A  common  assumption  m  using  this  approach  to 
model  soils  is  that  the  voids  are  either  spherical  or  flat  in  shape. 

O'Connell  and  Budiansky  (1974)  have  considered  the  effect  tnat  flat  cracks 
would  have  on  the  moduli  of  a  material  containing  such  voids.  The  equation 
they  obtained  for  the  bulk  modulus  of  such  a  material  is 
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bulk  modulus  of  material 

buxx  modulus  of  matrix  material 

Poisson's  ratio  of  material 

Poisson's  ratio  of  matrix  material 

crack  density 

crack  iength 

volume  of  material 
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The  Eqs.  (2.39)  were  aeveiopea  oy  considering  tne  cracks  to  contain  only  air. 
Tr.e  Eqs.  (2.39)  inaicate  tnat  a  sufficiently  iarge  crack  density  would  nave  a 
consiaeraoie  effect  on  the  material  properties,  wmie  tne  cracks  themselves 
may  be  of  negligible  volume,  as  the  pressure  is  increased  on  sucn  a  material, 
tne  cracks  would  close  and  tneir  effect  would  disappear. 

Other  researchers  have  considered  the  effects  of  spherical  voids  on 
material  benavior.  MacKenzie  (1950)  determined  tne  effective  bulk  modulus  for 
a  material  represented  by  a  matrix  containing  spnencal  voids.  The  term 
effective  refers  to  material  properties  which  are  descriptive  of  the  entire 
mass  being  considered.  The  porous  material  is  modeled  as  a  collection  o.f 
spneres  of  matrix  material,  each  containing  a  spnencal  void.  Under  this 
assumption,  tne  problem  reduces  to  that  of  determining  the  solution  for  one  of 
these  composite  spheres  with  a  uniform  radial  pressure  acting  on  its  boundary. 
The  term  composite  refers  to  the  material  composed  of  both  matrix  and  voids. 
The  expression  MacKenzie  obtained  for  the  effective  bulk  modulus  of  such  a 
material  is 
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(2.40) 


where  K  -  effective  bulk  modulus 

K  =*  bulk  modulus  of  matrix  material 
m 

G  -  snear  modulus  of  matrix  material 
m 

V  «  voiume  of  matrix  material 
m 

V  »  volume  of  composite  material 

The  Eq.  (2  40)  was  developed  under  the  assumption  that  air  is  contained  in  the 
voids.  Haanin(l970)  nas  determined  upper  and  low er  bounds  for  tne  effective 
bulk  and  shean  moduli  of  an  elastic  matrix  material  whicn  contains  spnericai 
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inclusions  ol  anocner  elastic  material.  Tne  upper  ana  lower  pounas  were 
determined  from  tne  tneorems  of  potential  energy  and  complementary  energy, 
ihe  bounds  determined  for  the  effective  bulk  modulus  coincided  and  this  result 
is  given  by 


(K  -  K  ) (4G  +  3K  ) c 
p  m  m  m 


4 G  +  3K  +  3  (K  -  K  )c 
o  p  m  p 


wnere  K  =  effective  bulk  modulus 


(2.41  ) 


*  bulk  modulus  of  tne  matrix  material 
K  *  buix  modulus  of  the  inclusions 
^ra  *  3iiear  modulus  of  tne  matrix  material 
c  -  the  column  fraction  of  the  inclusions 


The  bounds  Hashin  obtained  for  the  effective  shear  modulus  did  not  coincide. 
These  bounds  are  given  by 
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1  *  (1  '  G^/i  °  c 
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(2.42a) 


G-  [’  *  c2  - 


(2.42b) 


wnere 


lower  bound  on  effective  shear  modulus 


Gu  *  upper  bound  on  effective  shear  mobulus 
Gp  "  shear  modulus  of  tne  inclusion 
.ne  coefficients  y^  ana  y.^~^  are  determined  from  tne  equations 
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v  -  Poisson's  ratio  of  tne  matrix  material 
m 

v  -  Poisson's  ratio  or’  tne  inclusions 
? 

x  *  snear  stress  on  macroscopic  Sample 
Y  •  snear  strain  on  macroscopic  sample 


»  volume  concentration  of  spneres 

The  terms  x  ana  Y  appearing  in  Eqs.  (2.43)  anc  (2.44)  are  Known  from  tne 
boundary  conditions  usea  to  aetermine  the  limits  on  the  effective  shear 
modulus.  These  boundary  conditions  correspond  to  the  cases  wnen  the  surface 
tractions  ana  surface  displacements  are  known  on  a  sample  of  tne  composite 
material.  The  boundary  conditions  for  tne  case  wnen  surface  tractions  are 
known  are  given  by 


cr  ■  (x  ■  +  a,  y  ■  ♦  b) 
xy  ~ 


(2.40) 


The  boundary  conditions  for  the  case  wnen  tne  surface  displacement’s  are  known 


are  given  by 


u  ( x ,  y ) 
x 


(2.47a) 


u  (x,  y) 

/ 


(2.47b) 


The  oounds  determined  by  Hasnm  nave  been  successfully  used  to  approximate 
tne  effective  elastic  moduli  of  composite  materials. 

Some  spnencaj.  void  modexs  nave  been  developed  wmcn  account  for  tne 
pxastic  yiexding  of  tne  matrix  material.  Torre  U9ud)  developed  sucn  a  model 
ana  tne  resuit  ne  obtained  is  given  by 


P  -  t  xn( — ■ — 
m  a  -  i  ' 


,2.46a) 
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(2.48b) 


wnere  P  =*  pressure 


'i  »  yield  stress  of  matrix  material 
m 


V  =■  voiume  of  matrix  material 
m 


V  »  volume  of  composite  material 


A  problem  with  Eq.  (2.48a)  is  tnat  tne  matrix  material  is  considered  to  be 


fuily  piastic.  A  model  snouia  oe  aoie  to  aescnbe  elastic  as  well  as  plastic 


pnases,  wmcn  occur  for  botn  loaamg  and  unloading,  a  step  toward  including 


botn  elastic  and  plastic  pnases  is  to  prescnoe  a  work-hardening  rule  for  tne 


matrix  material.  ChaawicK  (1963)  developed  sucn  a  model.  Aitnough  this  model 


is  rigorously  derived,  certain  essential  parts  remain  in  integral  form  making 


it  difficult  to  use.  Carroi  and  Holt  (1972)  as  well  as  Chu  and  Hashin  (1971) 


took  an  approacn  wmcn  simplified  the  results.  Considering  the  same  spherical 


pore  geometry,  tney  derived  the  pressure-volume  relationship  for  the  composite 


material  by  temporarily  assuming  tnat  tne  matrix  material  is  incompressible. 


larroi  and' Holt  then  used  an  empirical  relationsnip  to  describe  tne  pressure- 


volume  relationship  for  the  matrix  material.  The  empirical  rexationsnip  used 


is  given  by 


1  s 

ri  ^  ^  >  '■*  =  ^ 


(2.49a) 


(2.49b) 


wnere  v  -  specific  voiume  of  composite  materia^ 


specific  voiume  of  matrix  material 


u  -  specific  internal  energy 

5 


Using  tne  metnod  outlined  aoove,  Carroi  and  Hoit  obtained  the  pressure-volume 
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a 


reiationsnip  of  a  mass  composea  of  an  iceaiiy  eiastic-piastic  matrix  material 
containing  voias.  This  is  given  by 
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(2.50(1) 


(2.50e) 


(2.50f) 


(2.50g) 


wnere  P  *  pressure 


G  =*  shear  modulus  of  matrix  material 
m 


Vq  -  initial  volume  of  composite  material 


V  -  volume  of  composite  material 


L  »  yield  stress  of  matrix  material 
m 


Tnere  are  two  problems  associatea  witn  using  Eqs.  (2.50)  to  represent  son 
oenavior.  First,  the  parameters  obtained  by  using  Eqs.  (2.49)  to  describe  the 
pressure-volume  relationsnip  of  tne  matrix  material  are  not  fundamentally 
related  to  tne  actual  behavior  of  tne  soil  grains.  Second,  soils  exhibit  a 
pronounced  reverse  yielding  during  unloading  wmch  is  not  predicted  by  Eqs. 
v2.ou;.  Bnatt  et  ai.  (1975)  attempted  to  remove  tnese  difficulties  by  making 
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tne  matrix  a  Monr-uoulomb  material.  The  yiela  criterion  for  the  matrix 
material  is  given  oy 


(1  +  D)  o1  -  a3  -  Y  -  0  (2.51 


where  a1  =  greatest  principal  stress 
=  least  principal  stress 
D  =  a  constant 

Y  *  yield  stress  for  tne  matrix  material 
.m 

The  results  obtained  by  Bhatt  et  ai,  using  tne  yield  criterion  given  by  Eq 
(2.51 ) ,  are  given  by 
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2  G  +  Y 


(2.52b) 


(2.52c; 


(2.52a; 


In  Eqs.  (2.52)  aix  terms  nave  tne  same  meaning  as  tnose  appearing  in 
Eqs.(2.5G).  The  parameter  x^  is  determined  from  tne  equation  given  oy 
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/men  tne  soil  is  aruoaded  in  tne  fuily  plastic  space,  tne  following 
relationships  nold 


Y  -^2.  Y  (2+D)  ^ 

p  *  -a  t(-a_) 3  -  i]  -  _ l  (cj 

D  Ot— 1  D  v  a' 

2  G  (ct-a*)  3  . 

+  -  m _  f  3+P  ,as  +  2, 

3  D (a-1  c  a 


(2.54a) 
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(2.54c) 


wnere  a*  =*  the  value  of  a  when  unloading  is  initiated. 

The  Eqs.  (2.54)  must  be  solved  numerically  to  obtain  tne  pressurerjvolume 
relationship  during  unloading.  This  is  done  by  choosing  values  of  c/a  >  1  , 
and  calculating  a  from  Eq.  (2.54b)  and  P  from  Eq.  (2.54a).  Some  problems  with 
tne  moael  just  described  have  been  recognized.  First,  the  predicted  high%- 
pressure  compressibility  is  often  too  low.  Second,  low  pressure  behavior  is 
not  adequately  represented.  Scnatz  et  al.  (1974),  modified  Bhatt's  model  to 
allow  for  the  curvature  of  tne  MonrmColoumb  failure  surface.  The  failure 
criterion  whicn  Scnatz,  et  al.  incorporated  into  Bhatt's  model  is  given  by 


Vo 
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(2.55) 


wnere  'i.  =*  yield  stress  for  -  0  condition 

0  max 

Y  »  ultimate  strength 

The  term  which  appears  in  Eqs.  (2.55)  is  descriptive  of  tne  rate  of. 

transition  from  low  to  high  yield  strengths.  Another  modification  wnicn 


Il 


5: 


Scnatz  et  al.  incorporated  into  Bhatt's  model  included  the  effects  of  flat 
cracKS  on  tne  bulk  modulus  of  tne  matrix  material.  Tms  modification  of  tne 
oulk  modulus  is  given  by 


K  ? 

—  =  1  -  Y  [1  -  — ],  <  Pcl  J 

m  cl 


(2.56a) 


K  "  Km  *  ^  -  ?d  *  (2’56b) 

wnere  Y  =  a  constant 

?Gl  =  pressure  required  for  complete  crack  closure 
K  »  effective  bulk  modulus  of  matrix  material 
«  buik  modulus  of  composite  material 
P  =  Pressure 

Tne  modification  given  in  Eqs.  (2.56)  nas  tne  effect  of  divining  tne  voids 
into  two  populations,  spnencai  voids  whicn  deform  according  to  Eq.  (2.55)  and 
flat  cracks  which  deform  according  to  Eqs.  (2.56).  These  modifications 
improve  the  preoictins  made  by  Bhatt's  model.  One  problem  with  Bhatt's  ana 
Scnatz' s  models  is  that  neither  allows  for  a  distribution  of  pore  sizes.  An 
approacn  to  account  for  the  pore  size  variation  in  actual  material  is  to  start 
with  the  ideally  plastic  spherical  pore  model  and  then  allow  for  each  sphere 
to  nave  a  different  porosity  with  the  requirement  that  the  total  porosity  is 
equai  to  that  of  tne  material  being  modeled.  Krener  and  Scnopt  (1*73)  have 
developed  such  a  model  which  considers  an  ideally  plastic  matrix  material. 
Tneir  resuxt  for  tne  pressure-volume  relationsnip  of  one  pore  is  given  by 
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(2.57b) 


where  P  31  pressure 


Ym  =*  yield  stress  of  matrix  material 
V  3  total  volume 
V!n  3  volume  of  matrix  maerial 
Vp  3  current  volume  of  pores 
Vp  3  initial  volume  of  pores. 

The  over.bars  in  Eqns.  (2.57).  denote  averages  taken  over  the  entire 
volume  of  material  being  considered.  The  pressure-volume  relationship 

for  the  entire  material  is  determined  by  evaluating  Eqns  (2.27)  for  all 

pores  present  in  the  material  under  consideration.  An  apparent  problem 

for  this  model  is  the  determination  of  the  pore  size  distribution. 

Other  spherical  void  deformation  models  have  been  developed,  but 

the  models  described  in  this  section  are  representative  of  work  which 

has  been  done  to  date  in  this  area. 


2. 3  Conclusions 

The  3oil  models  reviewed  in  the  previous  sections,  with  few  excep¬ 
tions,  have  only  considered  the  pores  within  the  soil  mass  to  contain 
air.  The  modification  of  many  of  these  models  to  represent  saturated 
3oil  conditions  is  straightforward  through  the  effective  stress  prin¬ 
ciple.  However,  many  situations  exist  when  the  soil  is  partly  saturat¬ 
ed.  The  degree  of  saturation  is  a  soil  parameter  usen  fo  describe  the 
amount  of  liquid  present  within  the  pores  of  the  soil  mass.  By  defini¬ 
tion  the  degree  of  saturation  is  given  by 
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S  =  (0  <  s  <  1) . (2.58) 

? 

where  3  -  degree  of  saturation 

Vp  *  volume  of  fluid  contained  in  the  soil  mass 
7p  =  vo  j.me  of  the  pores  contained  in  the  soil  mass. 

.he  definition  given  by  Eqn.  (2.58),  the  three  conditions  which  a 
constitutive  model  should  be  able  to  represent  are 

L.  S  ■»  0  (voids  completely  filled  with  air),  ...  .. 

2.  0  <  S  <  1.0  (voids  filled  with  an  air-fluid  mixture),  and 

3.  S  *  1.0  (voids  completely  filled  with  fluid). 

The  second  condition  presents  problems,  due  to  the  complexity  of  having 
an  air-water  mixture  present  in  the  pores.  One  problem  is  that,  a3  the 
pressure  is  increased  some  of  the  air  will  be  driven  into  solution. 
3ecause  of  this  and  the  compressibility  of  the  air  phase,  it  is  diffi¬ 
cult  to  predict  the  pore  pressure  resulting  from  the  application  of  a 
load.  If  the  pore  pressures  could  be  predicted,  the  principle  of 
effective  stress  could  be  used  to  model  the  partly  saturated  system. 

Phenomenological  models  have  been  used  a  great  deal  to  model  soil 
behavior.  It  would  seem  that  empirical  models  obtained  from  curve-fit¬ 
ting  mechods  are  undesirable  for  use  as  a  constitutive  model  represent¬ 
ing  soil  behavior.  These  models  should  not  be  expected  to  yield  rea¬ 
sonable  results  when  used  to  represent  conditions  which  deviate  greatly 
from  Chose  by  which  Che  model  was  calibrated.  They  also  provide  no 
understanding  as  to  the  actual  deformation  mechanisms  acting  within  the 
soil  mass.  Elastic  models  are  poor  representations  of  soil  behavior 
primarily  due  to  their  inability  Co  predict  unloading  behavior. 
Elastic-plastic  models  have  been  used  a  great  deal  and  provide  reason- 
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able  results  for  many  situations.  While  these  models  may  work  well, 
they  often  require  a  great  many  parameters  and  may  be  difficult  to  use 
in  practice.  Little  work  has  been  done  using  viscoelastic  models  for 
soils . 

Micromechanical  models  attempt  to  derive  constitutive  laws  from 
observing  the  actual  mechanisms  causing  deformation  of  the  microstruc¬ 
ture.  Thus,  as  these  deformation  mechanisms  are  more  fully  understood, 
a  better  understand  of  the  complex  behavior  of  soils  can  be  achieved. 
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CHAPTER  3 


ASSEMBLAGE  jF  SPHERES  IN  CONTACT 


1  Contact  Forces 


In  an  assemblage  of  spheres  whicn  are  in  contact,  tnere  are  tnree  types  of 


forces  whicn  may  act  on  tne  area  of  contact.  Two  of  these  are  forces  which 


act  m  directions  normal  and  tangential  to  tne  contact  area.  The  third  force 


is  a  torsional  couple  acting  on  tne  contact  area.  The  forces  are  shown  acting 


on  an  area  of  contact  in  Fijj.  3.1.  The  area  of  contact  results  from  the 


compression  of  one  sphere  upon  another.  The  remainder  of  this  chapter  will  be 


concerned  witn  the  normal  forces  whicn  act  on  the  area  of  contact. 


1  “I  Hertz  Solution  for  the  Pressure  Between  Two  Spherical  Bodies  in  Contact 


The  solution  for  the  pressure  between  two  spherical  bodies  in  contact  was 


first  determined  by  Hertz.  Discussion  of  this  solution  is  given  by  Timoshenko 


and  Goodier  (1951).  The  Hertz  Solution  will  be  reviewed  in  the  remainder  of 


this  section. 


Two  spherical  bodies  in  contact  are  snown  in  Fig.  3.2.  Here  R1  and  R^  are 


:ne  radii  of  spheres  1  and  2,  respectively.  Sphere  1  has  material  constants 


E  and  V  .  while  spnere  2  nas  material  constants  E_  and  V_.  The  x,  y  plane  is 


tangent  to  tne  point  of  contact.  The  z^  and  z^  coordinate  directions  are 


considered  positive  wnen  directed  from  the  origin  of  the  x,  y  plane  to  the 


centers  of  spheres  1  and  2,  respectively.  When  there  is  no  pressure  between 


-ne  bodies  tne  coordinate  directions  e  amd  are  given  by 


Z1  *  Rl“ [ r^] 


(3.1a) 


r  2  2 
z2  *  R2“fR2  “  r  1 


(3.1b) 
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Approximations  of  z1  and  z^  may  be  obtainea  by  performing  Taylor  series 
expansions  of  Eqs.  (3.1)  about  r  -  0.  These  approximations  are  given  by 


(3.^a) 

(3.2b) 


The  use  of  Eqs.  (3.2)  is  limited  to  cases  in  which  tne  distance  r  is  small  in 
comparison  to  R^  and  R^.  Addition  of  Eqs.  (3.2)  yields  the  following  equation 

r2  (R  -  R2) 

Z1  +  Z2  2  Ri  'R2  (3’3) 


As  shown  in  Fig.  3.1,  Eqs.  (3.3)  represents  the  distance  between  points  on  the 
surface  of  spheres  1  and  2  for  a  particular  value  of  r. 

If  tne  two  spherical  bodies  shown  in  Fig.  3-2  are  subjected  to  a 
compressive  force  F  directed  along  the  z1  and  z^  axes  such  that  tnere  is  force 
equilibrium,  the  bodies  will  make  contact  over  a  small  circular  surface.  The 
projection  of  tnis  surface  on  the  x,  y  plane  is  termed  the  region  of  contact. 

The  displacements  in  a  direction  normal  to  the  x,  y  plane,  of  points  lying 
on  the  surfaces  of  spheres  1  and  2  will  be  denoted  by  w^  and  w  ,  respectively. 
As  tne  spneres  are  pressed  together,  the  distance  between  two  such  points  on 
tne  region  of  contact  will  diminish  by 
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From  Eqs.  o .  5 )  ana  o.-tj  tne  following  relationsnip  is  ootamea  for  dj.i 
points  lying  on  one  region  of  contact. 


a  -  (w^  +  w2) 


r2^  +  R2) 

Z1  "  22  *  —  ^  R2 


(3.5) 


It  is  assumed  tnat  tne  radius  of  tne  region  of  contact  is  very  small  in 
comparison  to  tne  raun  ano  R^.  Then  wnen  considering  the  local 
deformation  witnin  the  region  of  contact,  the  spnere  may  be  considered  to  be 
repre^f|^ted  oy  a  naif-space.  This  enables  one  to  see  the  solution  for  a  point 
load  acting  on  an  elastic  naif-space  to  determine  the  displacements,  w.  and 
*2 >  within  the  region  of  contact.  The  geometry  for  the  problem  of  a  point 
load  acting  on  an  elastic  half-space  is  shown  in  Fig.  3.3.  The  solution  for 
the  dispiaceAat  z  «  0,  in  the  z-coordmate  direction,  is  given  by 


_  F(1-u2) 
z  -  0  ”  u  E  r 


(3.6) 


where  F  »  the  magnitude  of  the  point  load 
E  -  modulus  of  elasticity 
v  »  Poisson's  ratio 

The  Eq.  (3.6)  may  be  used  to  determine  the  displacement,  w,  on  the  plane  z  » 

0,  wnen  a  pressure  P(r)  is  applied  over  a  circular  area  on  this  plane.  For  a 
point  A  lying  within  a  circular  area,  as  shown  in  Fig.  3.4,  the  deflection,  w, 
at  point  A  may  be  determined  from  Eqs.  (3.6)  by  making  the  following 
substitution 


F  *  ff  P(s)  3  d<Cds  (3.7) 

B 

‘n  Eq.  (3.7),  the  integration  is  taken  over  tne  loaG  area,  B.  The 
displacement  of  point  A,  snown  in  Fig.  3.4,  is  given  by 
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s» 


I 


'J 


U-v2) 


2  =  0 


te. 


*  j  P  (  S  ,  p  )  dSQp 


(8.3) 


„  ,  7T  TT 

lH  tq.  1^3.8) ,  tne  angle  ip  ranges  from  —  to  j 


By  substituting  Eq.  (3.3) 


into  Eq.  (3.5),  tne  foilowiong  equation  is  obtamea. 


K,  .  K2 


H 


r/  P(s  ,<p)  asa<p 


(R1  *  ‘V  2 


2  R.  R  t  r 


(3.9) 


wnere 


Kl 


(l-vL  ) 


d-v2  ) 


Rc  -  the  area  contained  witnin  tne  region  ox'  contact 


The  pressure  distribution,  P(s),  is  chosen  sucn  that  Eq.  (3.9)  is  satisfied. 
The  pressure  distribution  which  satisfies  Eq.  (3.9)  is  an  exxiptical  cap  over 
tne  region  of  contact.  A  cross-section  of  tms  pressure  distribution  along 
the  chord  3C  shown  in  Fig.  3-9  is  shown  in  Fig.  3.5.  The  maximum  pressure, 


PQ(p),  axong  the  chord  3C  is  given  by 


1/2 


-r 


la  -  r  sin  <p  ;  ,  v-  \  *  p  1  \  j 


(3.10) 


wnere  P  =*  the  pressure  acting  on  the  center  of  tne  contact  region 


a  -  the  radius  of  tne  region  of  contact 
The  distribution  of  tne  pressure  axong  tne  cnora  5C  m  Fig.  3.5  is  given  by 
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P(S,4>) 


?o  (?)  2  2 

a  ""cos- 6  la  005  0 


(s  -  r  cos  ? )  | 

(-  |  <  a  <  tt/2  j 


(3.11  ) 


one  relationship  between  the  angles  6  aba  <p  is  given  by 


a  sin  9  *  r  sin  ? 


(3.12) 


It  is  useful  to  express  tne  variable  s  in  terms  of  an  angle  3.  This 
relationship  is  given  by 

S  -  r  cos  ?  +  a  cos  d  cos  3,  (0  1  3  -  it)  (3.13) 


Substition  of  Eqs.  (3.10),  (3-11),  (3,12),  and  (3-13)  into  Eq.  (9)  results  in 
the  following  equation. 


(Rl  +  Rj) 

2  h  h 


$)  sin^ 


r2 


3  d6d4>  * 


(3.14) 


The  expression  obtained  by  performing  the  integration  indicated  in  Eq,  (3,14) 
is  given  by 


(K, 


V 


2 

a  n 


2 

r  it 


(ri  :  v 

2R1  Ra 


I  r 


(3.15) 


Tne  Eq.  (3.15)  snows  that  Eq.  (3.9)  is  satisfied  by  an  elliptical  pressure 
aistribution  acting  on  the  region  of  contact,  provided  that  the  radius  of  the 


*2 


>v<> 
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contact  region,  a,  and  tne  displacement,  a,  are  given  by 


7T  ^  R2  (Kx  +  K2)  Pc 


2  0^+  R2) 


(3.16a) 


a  (K1  +  K2)  Pq 


(3.16b) 


The  constant,  P  ,  is  determined  from  the  static  equilibrium  of  one  spnere. 

For  equilibrium  of  one  sphere,  the  integral  of  tne  pressure,  P(s),  over  the 
region  of  contact  must  equal  tne  force,  P,  pressing  tne  spneres  together.  The 
pressure  distribution,  P(s),  symmetric  with  respect  to  the  center  of  tne 
contact  region  and  is  given  by 


P(r )  -  -2.  <  a2  -  r2 


,  (r  £  a) 


(3.18) 


The  condition  for  tne  equilibrium  of  one  of  tne  spheres  snow n  in  Fig.  3. <2  is 
given  by 


// 


P(r)  dA  *  2TT 


/  f 


rdr  *  F 


(3.19) 


Integration  of  Eq.  (3.19)  yields  the  following  value  for  P 


o  ,  2 

2ira 


(3.2o; 


The  displacements,  w1  and  w„ ,  normal  to  tne  region  of  contact  and  lying  within 
tne  region  of  contact  are  determined  from  Eqs.  (3.5).  These  displacements  are 
given  by 


TT  ^  P 


-  [2a2  -  r2]. 


(r  £  a) 


(3.21a) 
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b  ;v 


*  K  P  .  , 

2  o  . _  2  2. 

—  f2a  '  r  1» 


(r  <  a) 


(3-21 d) 


t 

Z\ 

ir 


s 

v’ 
v 
*  *. 


The  solution  obtained  for  the  two  spheres  in  contact  provides  information 
aoout  the  displacements  ana  stresses  occurring  on  the  surface  of  tne  spnere, 
witnin  the  region  of  contact.  It  does  not  provide  a  solution  describing  the 
stress  ana  displacement  fields  in  the  interior  of  the  sphere. 

3 . 3  General  Solution  to  tne  Axisymmetnc  Field  Problem  of  Elasticity  for  a 
Region  Bounded  by  a  Sphere 

In  this  section  tne  solution  of  tne  elastic  fiela  equations  for  a  spnere 
subject  to  eitner  axisymmetric  surface  displacements  or  surface  tractions  is 
obtained.  A  spnere  is  snown  in  Fig.  3.6  relative  to  both  tne  rectangular 
coordinates,  (x,y,z)  and  tne  spherical  coordinates  (p,p,9). 

-ihe  following  restrictions  are  imposed  on  the  spnere 

1 .  Surface  displacements  or  tractions  are  axisymmetric  with  respect  to 
tne  z~axis. 

2.  The  spnere  is  in  static  equilibrium. 

3.  Body  forces  are  negligible. 

The  approacn  taken  to  obtain  a  solution  to  this  problem  is  to  use  Boussinesq's 
solution  m  the  harmonic  function. 


£ 


u 


Tne  general  solution  for  a  region  with  torsionnfree  rotational  symmetry, 
in  tne  absence  of  oody  forces,  may  be  obtained  as  tne  sum  of  tne  two 
displacements  fields  given  by 


2G  m  -  7  $  (r,z) 

2G  “2  *  7  I*  1  (r,*)]  -  4  (I-v)  H  ( r,*)i 


(3.22) 

(3.23) 

wnere  UI  > ^  ”  displacement  vectors  in  rectangular  coordinates 
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7  *  gradient  operator 

e  ,e  ,e  =  amt  vector  in  tne  x,y,  ana  z  coordinates  directions 
x  y  z 

respectively 

2  2  2  2 
r  -  x  +  y  +  z 

G  »  snear  modulus 


v  =  Poisson's  ratio 

In  Eqs.  (3.22)  ana  (3.23)  $(r,z)  ana  v(r,z)  are  arbitrary  narmonic  functions. 
Henceforth,  the  solutions  given  by  Eqs.  (3.22)  and  (3*23)  will  be  referred  to 
as  tne  first  and  second  Soussinesq  solutions,  respectively . 

The  region  of  interest  is  bounded  by  a  sphere;  therefore,  it  will  be 
usefux  to  employ  spherical  coordinates.  Spherical  coordinates  are  related  to 
rectangular  coordinates  through  the  mapping 


x  »  p  sin<p  cos6 


y  -  p  sin<p  sine 


z  -  p  cos<? 


(3.24a) 


(3.24b) 


(3.24c) 


The  relationship  between  tne  displacement  components  in  spherical  coordinates 
ana  in  rectangular  coordinates  is  given  by 


u  »  3inq>  cose  u  ♦  smo  sine  u  *  cos<?  u 
s  x  y  z 


(3.25a) 


u,«  -  cos<p  cose  u  *  coso  sine  u  -  sm<p  u 
w  x  y  i 


(3.25b) 


u  -  sine  u  *  cose  u 
e  x  y 


(.3.250) 
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wnere  u  ,  u  ,  u„  »  components  of  tne  displacement  vector  m  tne  p,  ? ,  ana 
P  P  o 


3  coorainate  directions  respectively. 


u^,  u  ,  uz  -  components  of  the  displacement  vector  in  the  x,  y,  and 


z  coordinate  directions  respectively. 


From  Sqs.  (3 .2.2)  and  (3.25),  tne  displacement  components  m  spherical 


coordinates  for  the  first  Boussinesq  solution  are  given  by 


34  34  34 

2G  u„  *  sin4  coso  -r—  +  sm<p  9in0  t—  +  cosq 

0  ox  ov  a z 


(3.26a) 


34  34  34 

2G  ■  cos4  cosQ  -e—  +  cos4  sinQ  t"  -  sin4  t—  (3.26b) 

4  dX  oy  oz 


(3.26c) 


2G  u0  *  0 


The  combination  of  Eqs.  (3.23)  and  (3.25)  gives  the  following  displacement 


components  for  the  second  Boussmesq  solution. 


34  34 

2G  u,  *  sia4  cos9  ♦  sin4  sinQ 


*  C08<4 


il 

9z 


(3.27a) 


(3-4V)'f 


34  _  3f 

2G  *  cos0  cos®  +  cos 4  sinQ 


-sia9  (3-4v)4f 


(327. b) 


2G  »  o 


(3.27c) 


4  *  4(0,4) 


(3.26a) 


'?  -  f(p,ip) 


(3.28b) 


mmm 
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From  Eq.  (5.2baj  tne  partial  derivatives  of  tne  narmonic  function  <p  appearing 
in  Eqs.  ii. 26)  are  given  by 


ii  ii  3o  3i  li 

3x  a  3p  3x  3<|>  3x 

ii  ,  ii  i£  +  ii  ii 

3y  3p  3y  34>  3y 

34*  3<$  3p  3$  3<£ 

3z  3p  ~tz  +  *5$  3z 


(3.29a) 


;3.29b) 


(3.29c) 


From  Eqs.  (3.28b)  tne  partial  derivatives  of  tne  narmonic  function  appearing 
m  Eqs.  (3.27)  are  given  by 


3o  3p 

3p  3x 

3ip  3 p 

3p  3y 

3  Ip  up 
3p  3z 


ay  3tp 

dtp  3x 

3i p  a  $ 
3$  3y 

d  V  dtp 

dtp  3z 


(3.30a) 


(3.30b) 


(3.30c) 


Tne  gradients  of  p  and  <p  are  obtainable  from  Eqs.  (3.24).  These  gradients  are 


-*  3p  dp  dP 

vp  dx  x  dy  ey  3z 


sintp  cosd  e  *  sm<p  sin©  e  *  cost?  e 
x  y  z 


(3.31a) 


3<P  dtp  dtp 

<  m  — —  g  +  t—  g  +  p 

dx  x  dy  y  dZ  z 


cos<p  cos©  cosip  sin© 

p  ex  p  ey 


-  - — 21  g  (3.31b) 

Combining  Eqs.  (3.26).  (3.29),  and  (3.31)  gives  tne  following  displacement 
components  for  the  first  Boussinesq  solution 


(3.32a) 


sintp  o$ 

p  3(rnsP) 


(3.32b) 


71 


I 

%  , 
& 

8 

& 

fc: 

c 

tt 

1 

jy 

V, 

F; 


r* 


Coiaoirung  Eqs.  (3.27),  (3.30)  ana  (3.31  )  gives  tne  following  displacement 
components  for  tne  second  3oussinesq  solution 


2 G  u0 


COS<P 


(3-4^)f 


2 G  “-sin1? 


C03<P  T(^oi$7  *  ^-4V)f 


(3.33a) 

(3.33b) 


In  Eqs.  (3.32)  and  (3*33),  cos^  nas  been  cnosen  as  an  independent  variable. 

The  harmonic  functions,  $  and  4 1,  wili  oe  represented,  in  part,  by  spnerical 
.narmonics  wmcn  are  functions  of  coso. 

The  strainrdispiacement  relationships  referred  to  the  spnerical  coordinate 
system,  for  the  case  of  rotational  symmetry  about  the  z-axis  are  given  by 


c  * 

PC  3o 


£ 

PP 

1 

«  — 

P 

3<p 

♦  ^ 

P 

£ae 

+ 

0 

coc<p 

p 

u<& 

e 

Pd 

1 

2 

1 

L  P 

dUp 

3c?  * 

3u<? 

3p 

£  „ 
<po 

-  0 

£ 

dp 

-  0 

(3.34a) 

(3.34b) 

(3.34c) 

( 3 » 34  d ) 
(3.34e) 

( 3 - 34f ) 


From  tnese  strain-displacement  relations  and  the  displacement  components  given 
oy  Eqs.  (3.32)  and  the  strain  field  for  tne  first  Boussmesq  solution  is 


2G  £ 


OP 


7r  -  „  3in2t> 

2G  2 

0 


2 

3  ?  cosfl  3$ 

2  2 

3(cos<|>)  p  3(cos4>) 


.  1  3$ 
0  3o 


(3.35a) 

(3.35b) 
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« 


* 

tr‘ 


V. 

V. 


£ 


i 


$ 


2G  £ 


3®  cos®  o$ 


99  0  3p  2  3(cos®) 

p 


2G  £ 


sin®  3®  sin®  3  ® 

P®  *  2  3(cos®)  p  3p3(cos®) 

P 


e  -  G 
®0 


c  -  a 

®p 


(3.35c, 


(3.35a) 


( 3  -  35  e) 


(3.35f) 


Substituting  tne  aispiacement  components  given  by  Eqs.  (3.33)  into  the  strain- 
displacement  relations  given  by  Eqs.  (3*30  yields  the  following  strain  field 
for  tne  second  Boussinesq  solution 


2G  e 


PP 


2G  e 


®® 


32®  a® 

cos  ®  L  p  — -  n  2 ( 1  “2v)  J 
3p 

.  2  2  .2 

3w  cos  ®  ®®  sin  ®cos®  a  ® 

3P  P  3( cos® )  p  ^(cos®)2 


(3.36a) 


-  2 ( 1 ~2v) 


2  „ 

sin® a® 


p  3 (cos®) 


2G  eQ  -  cos® 


a®  cos  ®  3® 

3p  "  p  3(cos®) 


(3.36b) 

(3-360 


2G  cp<p  -  2(1-v)  -nf-°SV  -(^T®)  +  (1"2v)  3in<P 


ap 


-  3 in®  cos® 


^  3  2® 
3o3 (cos®) 


£  „  -  0 
®0 


e,  -  0 
0P 


(3.36a) 


( 3  •  36e) 


:3.3of) 


The  stress  fields  corresponaing  to  the  two  solutions  are  aeterminea  from  the 
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constitutive  equations  for  an  isotropic,  linearly  elastic  material.  Tnese 
constitutive  equations  are 


2G  e  +  Ae 
PP 


(3.37a) 


2G  e  +  Xe 
P<P 


(3.37b) 


2G  +  As 
a9 


(3.37c) 


■r 


(3.37a) 


(3  •  37e) 


2G  e  . 
p  o 


( 3  -  37  T ) 


w'nere  e  -  volumetric  strain  -  e  +  e  ^  +  eac 

PP  pp  09 


"(1-2v) 

Suostitution  of  tne  strain  field  given  by  Eqs.  (3-35)  into  the  constitutive 
relations  given  by  Eqs.  (3-37)  gives  tne  following  stress  field  for  tne  first 
Boussinesq  solution 


i3.36a; 


2  .2.  , . 

sin  p  o  p  cosp  dp 

2  ,  ,2  “  2  d(cosp) 

p  o(cosp)  p 


L  ££ 

p  ap 


(3.33b) 


1  3$  cosOi  3$ 

D  3p  2  3(cos<5) 


(0.33c) 
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S 1  ny  1  oy 

Pp  p  *•  p  3(cosy) 


,2 

a  £ 

dpo(COSy)  J 


0-3bd) 


o  »  0 
y  0 


(3-38e) 


a  -  0 

0p 


(3.38f) 


.Substitution  of  tne  displacement  components  given  Py  Eqs.  (3-36)  into  the 
constitutive  relations  given  by  Eqs.  (3*37)  yields  tne  following  stress  field 
for  tne  second  3oussmesq  solution. 


,2 

dp  op 

p  COSy  — r—  -  2  (1-v)  cosy  -r~ 
~  a  ap 

op 


sin f  dy 
p  3 (cosy) 


(3.39a) 


2  .2 

sin  ycosy  o  y 

p  3 (cosy) 


2  ♦  (1-2v)  cosy 


(3.39b) 


_  .  ay  i ,  2  2 

-2v)  cosy  ~~i  cos  ?  *  2Vsin  t>  J 


oy 

3( cosy ) 


(3.390 


, ,  „  ,  ay  ,  cosy  oy 

siny  i  (1“2v)  - —  +  2(1-0  -  T7 - r 

L  op  ‘  p  3 ( cosy ) 


.  a  y  i 

COS^3p3(cosO  J 


(3.39a) 


a  ■  0 
y  o 


( 3 • 39e) 


a  *  0 
op 


( 3  •  39  f ) 


)ne  term  of  tne  narmomc  functions,  y  and  y,  will  oe  determined.  These 
"unctions  satisfy  Laplace's  equation.  When  written  in  spherical  coordinates. 


IrV.1**** 


a 
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CM 15  5QudC10n  15 


1  3“3 


323  2  30  cot?  30  _L_  3^3 _ 

rT  P  oo  2  *P  2  ,  2  2  .  2  -02 

3p  H  p  p  3?  p  sin  P  39 


(d.uo; 


t: 

v  > 


A  particular  solution  to  Eq.  (3.40)  is  of  cne  form 


pa  ¥ 


(3.41  ) 


jr 

v 


wnere 


If  (cost®) 
a 


,V 


substituting  the  solucion  given  by  Eq.  (3-A1)  into  Eq.  (3-40)  yields  the 
following 

*2*  3* 

ain2T  2  co.»  srasay  .  -(n-l)  »  -  0  (3.42) 

The  Eq.  (3.42)  is  Legenare's  equation,  for  which  there  are  two  solutions. 
These  solutions  are  given  by 

¥•  =»  P  (cos$)  ,  (~  ao<  n  <  00 )  (3.42a) 

n  n  — 

-9q  ■  Qn  (cos?)  ,  (-  ®  <  a  <  oo)  (3.U2b) 


c»: 


The  solution  P  is  caiiea  the  Legendre  polynominal  of  degree  n. 
n 

m  is  caiiea  tne  Legendre  function  of  tne  second  kind  of  degree 
n 

solution  Q  contains  a  iogancnmic  singularity  cosip  -  ±  to  tne 
consideration.  Tne  solution  to  be  employee  is  P  .  An  equation 


The  solution 
n.  The 
prooiem  under 
defining  ?n  is 


M 

v 
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Pn  cos  ?) 


2nn!'  d(cos<t)a 


(cos^  $  -  x )n 


(3.^3) 


hencefortn,  tne  argument  or’  Pn  is  eosp.  Recurrence  rexations  for  tne  Leger.are 
polynomials,  appxicaole  for  all  values  of  n  are. 


P  -  ?  . 

n  n-1 


.Inn)  cosy  ?n  -  (n+1  )  Pn+1  -  n  ?n_1 


am"}  ?’  -  n  P  .  -  n  eosq>  ? 

n  n-1  *  n 


(3.44a; 


(3.44b; 


(3.44c) 


wnere 


n  d(cosq) 


The  narmonic  functions,  p  ana  p,  will  oe  of  tne  form  given  by  Eq.  (3-41), 
wnere  7  is  given  by  Eqs.  (3.42a).  3y  virtue  of  tne  first  recurrence  relation 
given  by  Eqs.  (3.44),  tne  functions,  p  and  y  are  given  by 


9  »pn1p  ,  (-  ®  <  n  <  ®) 
n  n 


~n-l  „ 

vn  -  p  pn  >  <-  -  <  n  < 


(3.45a) 


(3.45b) 


In  tne  Eqs.  (3-45),  tne  functions  ?n  and  y  represent  component  solutions.  A 
solution  to  a  particular  problem  character izea  by  specific  boundary 
conditions,  win  be  determined  by  superposition  of  tne  component  solutions 
ana  y  .  The  component  solutions  to  tne  first  and  second  Boussmesq  solutions 
win  nenceforth  be  denotea  as  wA^j  ana  respectively.  Using  Eq.  (3.45a), 

tne  aispxacement  3train,  ana  stress  fiexas  for  solution  (A  ,  are  given  by 


m 
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V 


;v. 


a 


ep 


=  0 


(3.460) 


Using  tne  expression  for  9  shown  in  Eq.  (3. 45b),  the  displacement  strain,  and 
stress  fields  for  solution  [CnJ  are  given  by 


S 


I  5 

!  .. 

i  ' 

• 

,  .  s 

f 


2G 

u 

P 

-(n+4-;4v) 
Gtn+1  )pn+1 

L  (n+1 )  Pn+1  +  n  Pno1  J 

(3.47a) 

2G 

i  1  3 

-  sm<p 

L  (n-3+4v)  P^+1  +  (n+4-4v)  P^  j 

(3.47b) 

Li  «a 

9 

(2n+1 )pn+1 

U0  " 

0 

(3.47c) 

2G 

(n+1 ) 

L  (n+1 ) (n+4n4v)Pn+1  +  n(n+4-4v)Pn_1  J 

(3.47d) 

"pp 

(2n+1 )pn+2 

2G 

c* 

1 

2 

L  (2n+1)  P^  n  (n+1)(n  -n+1+4nv)  PQ+1 
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(2n+1 )pn+2 

nn(n+1 ) (n+4n4v)  P  I 

nnl 

(3.47e) 

2G 

e00 

•  t  ("♦')  V,  '  en  > 

(3.47f) 

2G 

s  imp 

L  (n2+2nri+2v)  p^+1 

""pq 

,  ,  ,  ,  n+2 

(2n+1  )p 

+  (n+1  )  (n+4r-4v)  P'  ,1 
nr  1  J 

(3-47s) 

9  0 

=*  0 

(3.47h) 

£ 

op 

=■  0 

(3.47i) 

rr 

9  _ 

(n+1)  r , 

(a2+5n+4-2v)  P  +  n(n+4-4v)  p  ,  1 

a+1  n-1 

( 3 • 47  j  ) 

GO 

(2n+l)pn+2  L 

- 1  r. 

...  2 

(2n+l)on+2  [  'Q+1"2V)  *  a(o+l) (u+4-4v) 

Pn+1 

-  (2n+l )  P' 
n 

■ 

(3.47k) 
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id  “  l  <n*1)(1-2v)  Pn+1  -  P<  j 

P 


(3-^71) 


_sin®  2 

^»=  - n+2  L  (n+1  )  (n+4-s4u)  p’  +  (n  +2n-  +2v)  P^+1  j  (3.47m) 

(2n+l )p  . 


a  =  0 
<p0 


a  =0 
dp 


(3-47n) 


(3.«7o) 


The  boundary  conditions  to  be  considered  correspond  to  the  two  cases  when 

either  the  surface  displacements  or  surface  tractions  are  specified  on  the 

sphere.  For  the  case  when  axisymmetnc  surface  displacements  are  specified, 

tne  known  displacements  components  will  be  u  and  u..  These  quantities  may  be 

P  9 

represented  by 


u  -  u  (R ,  cosp)  on  S 
P  P 


u  »  u  (R ,  cos® )  on  S 
9  9 


where  S  -  the  surface  of  the  sphere 
R  -  the  radius  of  the  sphere 


(3,48a) 


(3.48b) 


s 
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ror  tne  case  wnen  axisiwtric  surface  tractions  are  specif  lea,  tne  known  stress 
icmponents  will  oe  ana  a  .  These  quantities  may  be  representea  oy 


a  =■  a  (R,  cos<p)  on  S 
PP  PP 


a  =  a  (R ,  cos® )  on  S 
P<P  P8 


(3.49a) 


(3.49b) 


To  appiy  tne  bounaary  conditions  given  by  Eqns.  (3.48)  or  Eqns.  (3.49),  these 
ooundary  conditions  must  be  represented  as  series  expansions.  For  tne 
displacement  boundary  conditions  given  by  Eqns.  (3.48),  these  series  have  the 


u  (R,  cos<p)  -  y  E  P 

p  rlo  u  " 


u  (R,  cos?)  *  sin?  )  nu  p* 


(3.50a) 


(3.50b) 


The  Eqn.  (3.50a)  is  a  series  expansion  for  a  Legendre  polynominal.  Eqn. 
(3.50b)  is  the  series  equation  for  an  associated  Legendre  function  of  the 
first  Kina,  of  degree  n.  The  significance  of  Eqn.  (3.50b)  is  that  the  cos#  » 

±  1 ,  the  series  expansion  is  equal  to  zero,  which  is  required  by  the  condition 
of  symmetry.  The  coefficients  in  the  senes  are  given  by 


(3n+1 ) 

— 5 — —  /  u-  ( R ,  cosu)  Pn  (cos<p)  d  (cos<p) 

-1 


(n-0,  1,  2,  3,...) 


(3.51a) 


(2nM)(n-1)1  .... 

- .  2(  ]  +n)|  J  u^  ( R ,  cos(j> )  si  nip  P„  (cos<p) 

d  (cos<p) ,  (n-  1,  2,  3,...) 


(3.51b) 


'.Vi.*  tLi  ^ 


v 
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£ 


In  Eqns.  (3.51),  me  superscript:  u  is  usea  to  aenote  tnat  surface 
displacements  are  tne  specified  boundary  conditions.  Expansion  of  tne  surface 
tractions  given  by  Eqns.  (3.49)  is  done  in  tne  same  way  as  tne  surface 

displacements,  except  that  the  coefficients  appearing  in  Eqns.  (3.50)  and 

a  a 

(3-bl)  are  denoted  by  £n  and  R^.  The  superscript  o  is  used  to  denote  that  the 
surface  tractions  are  the  specified  boundary  conditions. 

The  Eqns.  (3.50)  show  that  the  expansions  of  the  boundary  conditions  are 
m  terms  of  Legendre  poiynominals,  or  their  derivatives,  of  one  degree  n.  In 
orcer  to  meet  these  boundary  conditions,  the  quantities  being  specified  for 
solutions  [A  j  ana  [CnJ  snould  also  be  in  terms  of  Legendre  poiynominals,  or 
their  derivatives,  of  one  degree  n.  Solution  [AnJ  satisfies  this  requirement, 
but  solution  rcnJ  does  not.  Therefore,  the  solution  C j  is  formed  from  a 
linear  combination  of  solutions  [Anj  and  C cn] .  The  component  solution  [BnJ 
is  given  by 


L! 


t  <  • 


i 


* 


i, 

s 

< 

\ 

* 


V, 

V, 

V. 


a  jS 

I  l 

7 


[3  ]  -  (2n+1)[C  J  -  (n+4-4v)  [A  J 
n  n  nri 


(3.52) 


From  Eqn.  (3.52),  the  displacement,  strain,  and  stress  fields  for  solution 

[3  j  are 
n 


2G  u 

P 


2G  u 

<P 


(n-H  )  ( n-*-4-4v) 

n+1  pn+1 

P 


sinj) 

n+1 

P 


( n~3+4v) 


u 


0 


0 


(3.53a) 

(3.53b) 

(3.53c) 
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(n+1;  \j)n 

—  u  €■  3  “  .  "  —  •  , 

pp  n+1  n+1 

P 


(3.53a) 


l  (nn3  +  ^  Fj!j  -  (n+1  ) 


(n  -n+1 +3nv)  Pn+1  j 


(3.53e) 


L  (2n+1  )  (n+1 )  P+  1  +  (n-3+4v)  j  (3-530 


sin  2 

~ n+2  (n  +2n-1+2v)  P^+1 
P 


(3.53s) 


e  „  -  0 

ip  0 


E,  -  0 
0P 


(3.53n) 


(3.53i) 


(n+1  )  (n  +5n+4r;2v) 


(3-53J) 


1  2 

■^2  L  (n*1)(n  -n+1~2v  P  n+1~  (n~3+4v)  P'n]  (3.53k) 


—j  L  ( n+1  )  (2n+1  )  ( 1  ~2v)  Pn+1  +  (n~3+4v)  P^j  (3.53D 


sin<p  2 

~"n+2  (n  +2n~1 +2v)  P'+1 
P 


(3.53m) 


a  -  0 

<p0 


(3  - 5  3n ) 


o  -  0 
9  P 


(3.53o) 


The  solution  to  the  problem  under  consideration  will  be  aetermined  as  a 
suitable  combination  of  tne  component  solutions  [AnJ  and  [B  J.  The  form  of 
the  solution  is  taken  to  be 


00 
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[Si  -  z 

n- 

In  Eqn.  (3.54),  [S]  represents  tne  solution  to  tne  aisplacement ,  strain  and 

stress  fields.  The  coefficients  a  and  b  are  constants  of  superposition 

n  n 

wmcn  are  chosen  such  that  the  specified  boundary  conditions  are  satisfied. 

Some  values  of  the  constants  an  and  bn  may  be  determined  by  evaluating  the 
component  solutions  [A  ]  and  [BnJ  at  p  *  0.  These  solutions  become  singular 
at  the  origin  (p  =  0)  for  n  1  0.  The  solution  should  not  contain  these 
singularities,  as  the  condition  that  [S]  be  finite  at  the  origin  is  imposed. 
This  condition  requires  that 


bnCB* 


(3.54) 


an  ”  bn  “  °*  (n  ’  0) 


(3.55) 


Therefore,  the  solution  [S]  is  given  by 


L3  j 


L 

n-1 


[a  [  A  j 

-n  "n 


5_  [  B  ] 

•  .a 


(3.56) 


3-J.1  Solution  when  the  Displacements  are  Specified  on  the  Surface  of  the 
Spnere 

The  boundary  conditions  which  will  be  considered  are  the  displacement 

components,  u  and  u  ,  specified  on  the  surface  at  the  sphere.  Due  to  the 
P  ip 

restriction  of  axisymmetry,  these  boundary  conditions  can  be  represented  as 


up  -  Up  (S,  cos4,)  on  s  (3.57a) 

"  u$  (R>  C08<^  on  s  (3.57b) 
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wnere  3  -  tne  surface  of  the  spnere 


R  =  the  radius  of  tae  sphere 


Provided  tne  functions  u  (4,  cosp)  and  u^  (R,  cos<j>)  are  sufficiently  smooth, 


they  may  be  represented  by 


u  (R,  cos<p)  -  l  5  P 
H  n-0  1 


(3.58a) 


U(})  (R,  COS0)  *  sia<p  £  Up 

n-1  a  Q 


(3.58b) 


u  u 

The  coefficients  amd  are  given  by 


^2  '  •  ^  /  Up  (R,  cos<j>)  Pn  (cos<p)  d  (cosqj) 


(3.59a) 


(2n+1  )(nr1  )t 
2  ( n+ 1  )  1  • 


/'  u$  (R,  cos<p)  sin<p 


P^  (cos<?)  d  (cosqj) 


(3.59b) 


The  solution  to  the  interior  displacement,  strain  and  stress  fields  is  assumed 


to  be  given  by  Eqn.  (3.56).  Using  the  displacement  components  up  and  u^  of 


[S] ,  in  conjunction  with  Eqns.  (3.58)  yields  the  following  set  of  equations 


from  which  to  evaluate  a„  and  b  . 

n  n 


(1“2v)  R  b^  =»  G  5“ 


(3.60a) 


n  a  U  ,  -  (n+1  )  (nr;2+4v)  R^b  U  _ 
-<n-\  -nr-c 


-u  b1"0 

2G  4n  R.  . 


(n  -  1,  2,  3,...) 


(3.60b) 


*• 
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+  (n+5-^v)  R2b_^n_2  »  2G  n*  R1_!n 

^  n  -  1  ,  2 ,  3 .... ) 


(3.60c ) 


Solution  of  Eqns.  (3.60)  gives  tne  following  constants  of  superposition. 


G  L(n+5~4v)  £  +  (n  +  3n~2+4nv)  j 

_ n _ '  n 

~n=1  [3n+1.-r2V(2n+1)J  Rn+1 

(n  -  1  ,  2,  3,  • . . ) 

u  G  +  n  N^] 

"Q"2  [3n+l-2(2n+l)v]  Rn+1  ’  (°  *  °*  L>  2 . } 


(3.61a) 


(3.61b) 


The  displacement  fields  for  the  component  solution  [AnJ  and  [B  ]  vanish  for  n 
«  1 .  Therefore,  the  values  of  a^.  and  b^.  remain  undetermined.  The  solution 

ri  I  “I 

to  tne  problem  characterized  by  the  boundary  conditions  in  Eqns.  (3.57)  is 
given  by 


a  ,  [A  J 
^nnl  -n-1 


(3-62) 


In  Eqn.  (3.62),  the  constants  of  superposition  aU  .  and  bU  _  are  given  by 
Eqns.  (3-61).  The  coefficients  and  appearing  in  Eqns.  (3.61)  are 
evaluated  by  Eqns.  (3.59)  for  arbitrary  surface  displacements  of  the  form 
given  by  Eqns.  (3.57).  The  component  solutions  [A  ]  and  [B  _]  may  be 
determined  from  Eqns.  (3.^6)  and  Eqns.  (3.53)  >  respectively. 

3.3.2  Solution  when  the  Tractions  are  Specified  on  the  Surface  of  the 
Spnere 


The  boundary  conditions  to  be  considered  are  the  stresses,  a  and  0 

PP  ob 

specified  on  the  surface  of  the  sphere.  Because  of  the  restriction  that  the 
surface  tractions  be  axisymmetric,  these  stresses  may  be  represented  by 
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a  =  a  vrt.  coso )  on  3 
PP  PP 


(3.63a) 


o  =o  (R,  cosq )  on  3 
pip  P<? 


(3.63b) 


Provided  that,  the  functions  a  (R,  cosip)  and  0  (R,  coscp)  are  sufficently 

smooth,  they  may  be  represented  by 


app  (a«  CO80)  *  Z  C?  Pc 


(3.64a) 


ap(})  cos$)  »  simfr  £  p ' 

n-1 


(3.64b) 


a  a 

The  coefficients  £n  and  are  given  by 


f  Pn  + 1  ) 

— 2~~  $  °pp  cosip)  Pn  ( cosq> )  d  (cosip) 


(2n+1)(nn1 )1  .  . 

'  2 ( n-*- 1 T l"  '  ;  a  pcp(R'  cos*)  3in? 

P'  (cosip)  d  (cos<?) 


(3.65a) 


(3.65b) 


The  solution  to  the  interior  displacement,  strain,  and  stress  fields  are 


assumed  to  be  given  by  [S]  in  Eqn.  (3.56).  Using  the  stresses  a and  ap<p  of 

[S]  in  conjunction  with  Eqns.  (3.64)  yields  the  following  set  of  equations 

0  a 

from  which  to  evaluate  the  constants  an  and  bn> 


2  (1+v)  b_2  -  iQ 


n(nrl)  a_n_^  "  (n+1)  [ ( n+1 ) (n-2)~2vj  R  b_n_2 


(3.66a) 


R2  5°  , 
^n 


(n  -  1,  2,  3,...) 


(3.66b) 
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.(n-1)  a_n_L  +  (a2+2n-l+2v)  R2  b.n_2  -  R  'h.  ■ 


(n  »  1,  2 ,  3,...) 


(3.66c) 


The  Eqns.  (3.66)  are  not  compatible  for  n-1,  this  value  of  n  results  in  two 
equations  with  one  unknown.  3y  the  consideration  of  static  equilibrium, 
another  condition  may  be  obtained  which  will  render  Eqns.  (3.66)  compatible. 

A  sphere  with  symmetric  surface  tractions  about  the  z-axis  is  shown  in  Fig. 
3.6.  Under  symmetric  loading,  equilibrium  is  automatically  satisfied  in  the  z 
and  y  coordinate  directions.  Equilibrium  of  the  sphere  will  require  the 
following 


j  T  •  e  ds  -  0 
z 


wnere  T  -  stress  vector 


(3.67) 


e  -  the  unit  vector  in  the  z  coordinate  direction 
z 

In  Eqn.  (3.-o7)>  the  integral  is  taken  over  the  surface  of  the  sphere  S.  The 

stresses  wmch  will  be  known  on  the  surface  of  the  sphere  are  a  and  a  .  . 

PP  p4> 

Therefore,  tne  stress  vector  T  is  given  by 


i  -  a  e  +  a  e 
PP  P  PP  9 


(3.68) 


wnere  e  ,  e  -  tne  unit  vectors  in  the  p  and  <j>  coordinate  directions 
P  P 

respectively. 

Tne  unit  vectors  o  and  o  are  given  by 

P  p 


e  »  sin<p  cose  e  +  3in<p  sine  e  +  cosp  e 
p  x  y  z 

e  a  cosq  cose  e  +  cosp  sine  e  -  sin<p  e 
<P  x  y  T  z 


(3.69a) 

(3.69b) 


Substitution  of  Eqns.  (3.68)  and  (3.69)  into  Eqns.  (3.67)  gives  tne  following 
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condition  for  equilibrium. 


J  i  o  _ 

pp  cosg  -  a  s  mo  j  as  =  u 

o  v  J 


(3.70) 


The  stresses  a  and  a  may  be  expressed  m  tne  form  given  by  Eqns.  (3.64) 
PP  Pd 

The  surface  element  ds  is  given  by 


as  -  R  d0a(cos<p) 


(3.71  ) 


Using  Eqns.  (3.64)  and  (3.71)  in  Eqn.  (3.70)  yields  the  following  condition  of 
equilibrium. 


2  2ir  1  "  a  2 

R  /  /  l  cos<p  i  £  P'  (cos<p)  -  sm  ? 

0*1  n-0  ™  n 


l  n  P'  (cos<p)j  d0d(cos<p) 
n-1 


!3.72) 


Using  tne  recurrence  relation  given  by  Eqns.  (3.44c)  m  Eqn.  (3.70)  gives 


2  2ir  1  ®  °> 

R  /  / '  L  cos<p  l  E  P  (cosd)  -  i  n  n  P  (cos<p) 

0,  „  1 1  n»  #  i 

ml  n-0  n-1 


+  costp  i  q  n  Pn  (cos<p)  j  d0d(cos<p)  -  0 
n-1 


(3.73) 


Che  ortnogonality  relationship  applicable  to  Legendre's  poiynorainals  is 

1 


Pn(cO80)  PB(C03  0)  d(c0S<J)) 


m  -  n 


(3.74) 


2n+l>  n  *  a 


Noting  that  PQ(coS(p)  -  1  and  P^cos?)  -  cos$ ,  and  evaluating  Eqn.  (3-73)  in 
accordance  with  Eqn.  (3. 73)  gives  tne  following  equilibrium  condition. 


89 


S 


.v 


0 


— —  L  C,  -  £  q1  ,=0  (5.75) 

Tne  Eqns.  Q 3 . oo j  are  compatible  due  to  Eqn.  (3-74).  Solution  of  Eqns.  (3.66) 
yields  tne  following  constants  of  superposition. 

(n2+2n-l+2v)  f  +  (n+1)  (n2-n-2-2'j)  na 
,  J  »  n  ti 

®  “  ■  “  ■ '  “  j 

“a"1  2(n-l)  [n2+n+l-K2n+l)v]  ^ ~ 2 

(n-  2,3,4,...)  (3.76a) 

+  a  na 

b-n-2  - T3 - 2 - -  (n  »  0,  l,  2,...)  (3.76b) 

2  [n  -H*-l+(2n+l)v]  Ra 


The  coefficients  of  superposition  a_1 ,  a_ ana  b_1  remain  undetermined 
oecause  tne  stress  fields  vanisn  for  the  component  solutions  [A_1 j ,  [A_2],ana 
[3_^j.  The  solution  to  the  problem  characterized  by  the  boundary  conditions 
given  in  Eqns.  (3.63)  is 


C-‘ 

V’ 


s 


'r< 

i 


[S] 


a 

a-n-l 


fA-n-J 


+  Z 

n*0 


b 


a 

-n-2 


(3.77) 


m  Eqn.  (3.77)  the  constants  of  superpositions  are  given  by  Eqn.  (3.76).  The 
coefficients  ^  and  appearing  in  Eqns.  (3.76)  are  determined  from  Eqns. 

(5.65).  The  component  solutions  [A_  . ]  and  [B_  . J  may  be  determined  from 

~riwi  -snr-s  i 

Eqns.  (3.46)  and  Eqs.  (3.53).  The  solution  given  by  Eqns.  (3.77)  is  tne  same 
as  tnat  determined  by  Sternberg,  Eubank,  and  Baaowski). 

3 .4  Solution  to  Specific  Elasticity  Problems  for  a  Spnere  Subject  to 
Axisymmetric  Surface  Displacements  or  Surface  Tensions 
The  results  of  Sections  3.2  and  3-3  will  be  used  to  determine  elastic 
solutions  for  a  sphere,  subjected  to  a  number  of  particular  boundary 


oonai cions . 


Tne  oouncary  conditions  wnic.n  wi*l  oe  considered  for  tne  spnere 


are  as  follows: 

1 .  The  surface  displacements  resulting  from  tne  contacts  witn  two 
adjacent  spneres  along  an  axis  of  symmetry. 

2.  The  surface  tractions  resulting  from  tne  contacts  witn  two  adjacent 
spheres  along  an  axis  of  symmetry. 

3.  A  uniform  radial  pressure  applied  over  the  entire  surface  of  tne 
spnere. 

The  oounaary  conditions  listed  above  satisfy  the  restrictions  imposed  on  the 
general  solution  given  in  Section  3 • 3 •  Therefore,  the  results  of  Section  3*3 
may  be  used  to  determine  solutions  to  the  problems  characterized  by  the 
boundary  conditions  listed  above. 

3.4.1  Solution  when  tne  Surface  Displacements,  Resulting  from  Three  Sphere  in 
Contact  Along  an  Axis  of  Symmetry,  are  Known 
The  surface  displacements  on  the  region  of  contact  resulting  from  pressing 
two  spheres  together  are  provided  by  Hertz  contact  theory.  It  is  assumed  that 
tne  surface  displacements  on  tne  region  of  contact,  parallel  to  the  region  of 
contact,  are  negligible.  The  surface  displacements  perpendicular  to  the 
region  of  contact,  for  the  case  of  two  spheres  m  contact,  are  given  by  Eqns. 
(3.21  ). 

Tnree  spheres  in  contact  along  the  axis  of  symmetry  are  shown  in  Fig.  3-7. 
The  displacement,  u  ,  for  the  region  of  contact  is  assumed  to  be  given  by 
Hertz  contact  theory.  Therefore,  Eqns.  (3-21)  are  used  to  obtain  tne 
following  displacement  boundary  conditions  for  the  center  sphere  snown  in  Fig. 

3-7. 
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U_  (R,  C03<j)) 


3 (1-v  )F  2  .  ,  , . 

— ^ — ? -  [2a2  -  R  sin24>] , 

8  a  E 


(7T  -  4>'  <  <j>  <  tt) . (3.77b) 


where  F  -  total  force  transmitted  through  the  contact 


a  -  radius  of  the  region  of  contact 


<p'  -  angle  defining  the  region  of  contact 


R  »  radius  of  sphere 


E  -  elastic  modulus 


Using  Eqns.  (3.25),  the  spherical  displacement  components,  u  and  u .  .  are 

P  <P 


u  -  u  cos* 
P  z 


(3.78a) 


u  -  u  Sind) 
<pz 


(3.78b) 


The  boundary  conditions  given  by  Eqns.  (3-77)  and  the  development 


components  in  Eqns. (3.78)  can  be  combined  to  yield  the  following  boundary 


conditions  in  terms  of  the  spherical  displacement  components,  u  and  u  . 

p  <p . 


Mu-OF 


3  a  E 


[2a2  -  R23in2<J>]  cos $,  (0  <  j>  £0  *0 


up(R,cos<f>) 


(3.79a) 


2 (I— v  j F 

- 5 -  [2a2  -  R2sin2$]  cos$,(tt  -  <  $  <  tt) 

N  8  a  E  “  “ 


U^(R,COS0) 


“ — i  ^  [2a2  -  R23ia$]  sin4>,(0  <  <P  <<p'  ) 

8  a  E  “ 


[2a2  -  R^sia24>]  sia$,(7T  -  <p '  <  <p  it) 
3  aJ  E  ~ 


(3.79b) 
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To  determine  the  displacement,  strain  and  stress  fields  for  a  sphere  whose 
displacement  boundary  conditions  are  given  by  Eqns.  (3.79),  it  is  necessary  to 
evaluate  the  constants  and  given  by  Eqns.  (3.61).  In  this  case,  these 
constants  are  given  by 


-u  ^  3(2n+l)  (1-0 
16  a3  E 


I"  cosCtt-^1) 


L 


2  2 

[2a  -  R  sin<?>]  cos$  P  (cos<J>)  d(cos<{>) 

a 


1 

^cos<P ' 


2  2  2 

[2a  -  R  sin  4>]  costp  P  (cos<j>)  d(cos<P) 

a 


1 


(3.80a) 


.  COS  (TT-4>  *  ) 


u  _  3(2n+l)(l-v  )F 
a  16n(n+l)a3  E 


[-L 

*/'■ 


[2a2  -  R2sin2<J>]  sin2<(>  P  (cos(J>)d(cos<j>) 

n 


[2a2  -  R2sin2$]  sin2s>  P'  (cos<j>)d(cos<|)) 


cos<J> ' . (3.80b) 


The  expressions  for  the  constants  given  in  Eqns.  (3.80)  may  be  simplified 
since  the  functions  contained  in  the  integrals  are  odd  or  even  functions, 
depending  on  the  value  of  n.  These  functions  are  odd  for  an  odd  value  of  n, 
and  even  for  an  even  value  of  n.  Therefore,  the  constants  E^  and  n^,  will  be 
zero  for  odd  values  of  n.  In  view  of  this  the  constants  are  given  by 


ru  -3  (4n+l)  (1- 

52“'~T77 


u)F  f  ,  , 

-  /  [2a2-R2si 

*'COS0' 


in  4>]  cos 4*  P2n(cos4>) 


d(cos4>),  (a  -  0,1,2 . )  .  .  .  (3.81a) 


,u  ,  3  (4n+l)  (1-v 


Zn 


16n(2n+l)a 


^  f  [  2a2-R2sin24>^  *iaZ*  ?'2n(cos<t>) 

*COS<p’ 


d(cos<?) ,  (n  -  1,2,3 . )  .  .  .  (3.81b) 
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Evaluation  of  the  integrals  contained  in  Eqns.  (3.81)  yields  the 

•  ,  u  u 

following  expressions  for  ^2q  an<^  ^  . 


-u  -3(l-v2)F  (2a4-R2)sln^(j>’ 


8  a3  E 


-u  3(1- 


vi  d^-p  f  2  ( 2n+l . 

;  [2sin|, _  u  — 

8  z  E  (  1  ( 

(4n+l)P2n(cos4)  2nP2a_2(cos4)  * 


2  4 

1  (1-cos  d>'' 


(2n4*l)P2a+2  (cos<fr ' ) 


(4n+3) (4n-l) 


6P2n+4(c0S<n 


(4n-l) 


6  cos<j) *  P2n+3(cos<j>') 


(4n4-7)  (4n+5)  (4n+3)  (4o+5)  (4n+3) 


+  3[ 

+  £cc 


8  2^,1  P2n+2(cos<;),) 

(4n+7)(4n-l)“  C°S  *  (4n+3) 


3A ,  18  cos<fr 

+  COS  <p  -  />  r  i 

(4o+5) (4a- 


2a+l(cos0*) 


+  6  [  cos2$’  - 


(4n+5) (4a 


.  I  18cos<&'  3. 

+  (4a+3)  (Zn-5>~  c03  * 


_ 1  1 

in-3) J 

'  ]  P2n-1 


(4n+l)  P2a(cos0') 


(4n+3) (4n-l) 


(cOS<j>'  ) 


+  3  f  3  coa2cb'l  P2n-2(COS0,) 

3  (4n+3)(4n-5)  "  *  (4n-l) 


6  cos$'  P2a-3(cos^f)  6  P2a_^(cos<t 


(4n-l) (4n-3) 


—  I 

a-5)  )  ’ 


(3.82a) 


(4n-l) (4n-3) (4a 


(n  *  1,2,3 . )  ....  (3.82b) 
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3(l-v  )R  F  - (4n±l) —  f  [2sin2<J>'-ll  [P  (cos4>') 

8  a3  E  (2“+l)  C4^-!)  i  L  iJli2n-2UOSV  ; 

. -  2  P2n+2(C0S4>,)  .  2  “*♦'  P2n+l(cos£p,) 

-  P2a(coa®  ) ]  (4n+3) (4a+1)  +  (4n+l 


+  f  (4nf37(4^~  *  COs24>f  P2n(coS<{>,) 


4  (4n-l)cos<p'  P2n_^(cos0r) 


(4n+l) (4n-3) 


lCcos$  j  r  2  ,  6  1 

+  [  C°S  0  “  (4n+l)  (4n-5)  J 
2  cos#’  P-  ,(cos<j>')  2  P  (cos4>*)  ) 

p2n-2<“a*'>  + - wdf3 - *  -(Sayrggy  | 


(2n+l) 


,  (n  -  1,2,3 . ) . (3.82c) 


The  general  solution  [S]  given  In  Section  3.3  may  be  modified  so  that 
it  does  not  include  the  zero  terns  which  occur  for  odd  values  of  n. 
Then,  the  solution  to  the  problem  characterized  by  the  boundary  condi¬ 
tions  given  in  Eqns.  (3.79)  is 


!S1  •  £,  -2q-1  ‘W  +  L  b2n-2  tB-2n-2l 

n*l  n«o 


(3.33) 


The  constants  of  superposition,  a_2n_]_  b  2n-2’  are  S*ven 


-2n-l 


G((2n+5-4v)  Co  +2(2n2+3n-l+4nv)n“  ] 
4ii  zn. 

[6n+l  -  2(4n+l)  v] 


.  .  .  (3.84a) 


-2n-2  [6n+l-2(4n+l)v] 


(3.84b) 


.'7*;  . 


't\V 


T  .  '  .*  ■ 


The  component  solutions,  tA_2n-l^  aQ<*  ^®-2n-2^  ’  oa^  be  determined  fr°® 
Eqns.  (3.46)  and  (3.53> ,  respectively.  The  component  solution  tA_2n_]J 
is  given  by 


2G  »p  -  -  2np2n_1P2n . (3.85a) 

2 G  -  sin^p211-1?^ . (3.85b) 


U0*° . 

2G  £  -  2n(2n-l)p2n“2  P~  . 

PP  2n 

20  %  '  p2°"2tP2n-l  *  2n(2a-l)  P2nI 

2G  e00  *  “  P  P2a-1 . 

2G  ep<d  *  “  (2n”l)  si a<^p2n'2  ?2n  *  * 

£<t>@  ■  0  . 

£ep  *  0  . 

app  -  2n(2n-l)  P2n  1  ?ln  . 

°<p<p  *  P  ^P2n-1  ~  2n(2n'’1-)  P2n^ 


(3.85c) 

(3. 85d) 

(3. 85e) 

(3.85f) 

(3.85g) 

(3. 85h) 

(3.85i) 

(3.85j) 


(3.85k) 
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p2n-2  pi 

2n-l 


(2n-l)  p^n_2  3 in<J>  P' 

2n 


0 


(3.851) 

(3.85m) 

(3.85a) 

(3. 85o) 


The  compoaenc  solution  [B  ]  is  given  by 


2G  up  »  -  (2n+l)(2n-2+4v)  p2n+1  P2n . (3.86a) 

2G  .  (2n+5-4V)  p2a+1  9i-a<J>  . (3.86b) 

u9  "  0 . (3.86c) 

2G  £pp  -  -  (2n+l)2(2n-2+4V)  p2a  ? 2n .  (3.86d) 


2G  -  -  P2n£(2n+5-4V)  p^+1  -  (2n+l)  [  (2n+l) 2 

♦  2(n+l)(3-4v) ]  P2nj  .  .  .  (3. 86e) 
2G  £g9  -  -  P2n£(4a+3)(2n+l)  P2n  -  (2n+5-4V)  P'2n+1J  .  .  .  .  (3.86f) 
2G  £p$  "  f  4n(n+l)  “  l+2v|  p2n  sin$  P^Q . (3.86g) 


( 3 . 86h ) 


Se 

£9p 

a 

PP 

a09 

ap<j> 

a$6 


0 


0  .  (3. 861) 

-  (2n+l )  [ (2n+l)(2n-2 )  -  2v]  p2n  P2n . (3.86j) 


P  [(2n+l)(4n  +10n+7-2v)  P2n+(2n+5-4V)  P^n+1]  .  (3.86k) 

p2n[(2n+5-4v)  P£n+1  -  (4a+3)(2n+l)(l-2V)  P2n]  .  .  (3.861) 


(4n2+4n-l+2  )  P2*1  sin<P  P£n . (3.86m) 

0  .  (3.86n) 

0  .  (3. 86o) 


3.4.2  Solution  when  the  Surface  Tractions  Resulting  from  Three 
Spheres  in  Contact  Along  an  Axis  of  Symmetry  are  Known 
The  surface  tractions  resulting  from  the  contact  of  two  spheres 
are  known  from  Hertz  contact  theory.  On  the  contact  surface,  the 
normal  stress  is  known  from  Eqn.  (3.18).  For  the  problem  of  three 
spheres  in  contact,  as  shown  in  Fig.  3.7,  it  is  assumed  that  the  normal 
stress  on  the  contact  surfaces  are  given  by  Eqn.  (3.18).  Therefore, 
the  stress  boundary  condition  is  given  by 
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a  tR,  CosP) 

pp 


3RF  ,  2*  2aI  . 

- 5-  [cos  p  -  cos  p  ] 

2:ra 


1/2 


(0  <  p  <  $ '  ;  TT-p  <  <J>  <  TT) 


(3.87) 


where  F  ■  force  transmitted  through  the  contact 

a  *  radius  of  contact  region 
P’  *  angle  defining  the  region  of  contact 
R  *  radius  of  the  sphere 


In  order  to  use  the  solution  developed  in  Section  3.3,  the  stress 
components  on  the  surface  of  the  sphere,  O  and  a  ,  must  be 
determined.  These  stress  components  are  given  by 

2 

<3..  9  <J  cos  $ . (3.88a) 

pp  zz 

o  ■  a  sinp  cosp . (3.88b) 

pp  zz 

From  Eqns.  (3.87)  and  (3.88),  the  desired  boundary  conditions  are  given 
by 


“3RF  2  2  ^/^  2 

0^(8, cosp)  ■  - [cos  p  -  cos  p  ]  cos  p, 

2iTa 


PP 


(0  <  4>  <  4>'  ;  ir-0'  <  p  <  tt) . (3.89a) 


3RF 


1/2 


2  2 

cr  (R.cosp)  «  — — r  [cos  1  -  cos  p']  sinp  cosp  , 
P<P  2iraJ 
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To  determine  the  displacement,  strain  and  stress  fields  for  the  inter¬ 
ior  of  a  sphere  whose  surface  tractions  are  given  by  Eqns.  (3.89),  it 
is  necessary  to  determine  the  constants,  £a  and  Hn»  S^ven  ^y  Eqns. 
(3.65).  For  the  boundary  conditions  in  Eqns.  (3.89),  the  constants  are 


given  by 


cos  (m-— $ ' ) 


_a  3(2n+l)RF 
3 


4ir  a 


m  ”  T 

f,  2 

I  [cos  $  - 


,  1/2  2 

cos Vl  cos  <|>  Pn  ( cos 0)  d(cos<p) 


+  [cos2<j>  -  cosfy]  cos2$  Pa  (C0S4>)  d(cos^) 

*tos<|) ' 

(n  -  0,1,2, . . . . 


cos(  - <p ') 


2... 1/2  .  2, 


3(2n+l)KF  I [cos  9  -  cos  $' ]  sin  Icostp  P^cos^)  d(cos4>) 

4n(n+l)ira  L  J 

-1 


7  “ 


2  2  1/2  2 
cos  $  -  cos  <j>' ]  sin  tficosipP^  (cos<j>) 


d(cos <P)  (n  ■  1,2,3,...)  ....  .(3.90b) 


The  above  integrals  for  and  n°  may  be  simplified  since  the  functions 

contained  within  the  integrals  are  odd  or  even  functions,  depending  on  the 

value  of  n.  These  functions  are  odd  for  odd  values  of  n  and  otherwise  even. 

Therefore,  the  constants  and  n°  may  be  determined  as  follows 

n  n 


OTnjtTrw 


W.t'« 
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e 


9 

£ 

f 


& 


$ 

K 


I 


1 


£ 


I 


r  [co 

*/cos4>' 


1/2 

*  -3(2n+l)RF  I  [coa2^  .  cos^cf)1]  cos2(J>  Pn(cos®)d(cos<J>) , 
_n  2tt  aJ 


(n  -  1,2,3 . )  .  .  .  (3.91a) 


2n  3 

2n(n+l)7r  a 


r  [co 

-'costp' 


1/2 


a  _  3(2n+l)RF  /  r  2X  2A1,  ' '  .  2X  „  , 

/  [cos  <j>  -  cos  <j>  ]  sin  <p  cos<J>  P2n(cos<j>)d(cos4>) , 


(a  -  1,2,3, - )  .  .  .  (3.91b) 


Evaluation  of  the  integrals  in  Eqns.  (3.91)  yields  the  following 


_a  ,  a 

eqations  for  §2n  and  n2n. 


5? 


-3FR 
0  2-ira3 


3  7  4 

sin  (ft1  +  cos  (fr'sinfl'  _  cos  <p '  ^ 


( l-fsin^' j 
V  cos<j)  ’  J 


.(3.92a) 


f  ( 4n*H)  3FR  (-l)n(4n-2m) ! 

2n  2tt  a3  o-o  22nm!  (2a-m)  !  (2n-2o) ! 


-sin<t> 


(2n-2art-2)  icosj)' 


2  (m-nrf2) 


(2n-2mf4  22  (n"ni+:L)  [  (n-ofl) !  ]2  (2n-2a+4) 


2! 


■r, 

V. 


£ 


£ 


E 


n-m 


2  ,  ,  2  (n-m-' 


+  sinfl'  -  (2m-2mrf2) !  (k! )  cosd) 

(2n-2mH)  J-  2 («*+!), ,2 
k«o 


[  (n-nrt-1) !  ]  (2k+l) 


m-k+1) 
!k+l)!  . 


(n  —  1,2,3,....)  .  .  .  .  . 


. (3.92b 


.a  _  (2n-2) ' 3FR 

l-nia3 


2n 


l-i 


(-l)m(4n-2m) ! 


2m 


27r(2n-l)  !a  l  2  m! (2n-m) ! (2n-2m) ! 


sin£>' 


(2n-2nrt-l) 
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2  (a-nri-1) 


f 2a- 2m+2) !  cos'* ' 


l+sind 


22  (n-nri-1)  !  ] 2  (2n-2sri-l) 


COS<J>  1 


(2n-2a+2)!  (k!)  costt)1 


2 (n-m-k+1) 


(2n-2m+l)  ^  22  (a"m"k+1)  [  (n-nri-1)  !  ]2(2k+l) 


*1)  ‘ 

; 


+  (2n-l)  Z 


(-l)m(4n-2m-4) 


22(a~i)mi(2n-m-2)!(2n-2m-2)!  (2n-2m-l 


sin<J>* 


J4(n-m) 

_ 3in<t  * _  (2n-2nri-2)  (2n-2nri-2)  !  costfe ' _ 

(2n-2nrt-l)  (2n-2m-l)  cos2<|> '  22  (n"nH‘1)  [  (n-nri-1) !  I2(2n-2iiri-l)  (2n-2m-l) 


(l+sinift'N 
cos$'  J  ~ 


(2n-2nri*2)sind 


(2n-2m+l) (2n-2m-l)  ^ 


_  2  (n-m-k+1) 

(2n-2nri-2)  !  (k!)  costft1 
02  (n-m-k+1)  r  \  1 \ , 


[(a-nri-1) !]  (2k+l) 


3- 

i 


2n(2n-l) !  ^ 
(2n+l) !  ^2n  ’ 


(a  ■  1,2,3,....)  ....  (3.92c) 


The  general  solution  [S]  given  in  Section  3.3  may  be  modified  so  that 


it  does  not  include  the  zero  terms  which  occur  for  odd  values  of  n. 


Then  the  solution  to  the  problem  characterized  by  the  boundary  condi¬ 


tions  in  Eqns.  (3.89)  is  given  by 


Z  a-2n-l  ^-2n-l]  + 


’  b-2n-2  ^n^ 


.  (3.93a) 


(4n2+4n-l+2v)  +  2(2n-l) (2n2-n-l-v) 


-2n-2 


2 (2a-l)  [4n2+2n+l+(4n+l)v]  R2n'2 


.  (3.93b) 


(3.93c) 


i 


} 

'* 

£ 


I  ft 


a  ,  ^2n  +  2n  n2n 

2n"2  2[4n2+2n+l+(4n+l)v]R2ri 


The  component  solutions  [A  and  [B  2n  2^  are  8iv«Q  by  Eqns.  (3.85) 

and  Eqns.  (3.86),  respectively. 


in 


v 

I  jO 


3.4.3  Solution  £or  a  Sphere  Under  the  Action  of  a  Uniform 
Radial  Pressure 

A  sphere  subjected  to  a  uniform  radial  pressure  is  shown  in  Fig. 
3.8  .  The  boundary  conditions  for  this  sphere  are 


a  (R,cos<p)  »  -  a  ,(0  <  <p  <  IT) . (3.94a) 

DP  4  ~ 


I 


ap(p  (R, cos4> )  -  0  ,(0  <  4>  <  7T) . (3.94b) 


where  *  the  magnitude  of  the  uniform  radial  pressure 


To  determine  the  displacement,  strain  and  stress  fields  for  the  problem 


$ 

$ 


I 

•• 

.U' 


characterized  by  the  boundary  conditions  given  in  Eqns.  (3.94),  the 

_a  a 

constants  E  and  n  must  be  determined  from  Eqns.  (3.66).  For  these 
n  n 

boundary  conditions,  the  constants  ^  and  r)°  are  determined  as  follows. 

n  n 


.0  - ( 2n+l ) 


/;• 


Pn(cos4>)  d(cos<p)  , 


(a  -  0,1,2 . ) 


(3.95a) 
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‘Y.*V 


nj  -0,  (a  -  1,2,3 . ) 


(3.95b) 


EvaluaCing  Che  integral  in  Eqna.  (3.95)  yields  the  following  for  the 


constants  n  and  • 


£  u  ,  _  a 
^o  u 


a 

S  "  0  .  <n  -  1.2,3, - ) 


Hn  -  0  ,  (n  -  1,2,3, - ) 


(3.96a) 


(3.96b) 


(3.96c) 


In  Eqns. (3.96) , the  superscript  <JU  has  been  used  to  denote  that  the 

constants  were  determined  for  Che  case  of  a  uniformly  applied  pressure. 

a  a 

The  Eqns.  (3.96)  show  chat  Cq  and  Hq  are  equal  to  zero  for  n  ^  1. 

Due  Co  this,  Che  solution  [S]  reduces  to 


where 


b-2  [B-2] 


2(1+m) 


(3.97) 


The  component  solution  is  determined  from  Eqns.  (3.54).  Substi¬ 

tuting  solution  [ inCo  Eqo*  (3.97)  yields  the  following  displace¬ 
ment,  strain  and  stress  fields  corresponding  to  the  solution  [S]. 


(l-2v)  a  p 

_ u 

2G(l+v) 


(3.98a) 


106 


U/j)  *  0 .  (3.98b) 


ue  -  0 .  (3.98c) 


— (1— 2v )  a 

:pp  *  2G(1-K>) 


(3.98d) 


-(1-2V)  cr^ 

2G(l+v) 


(3. 98e) 


-(l-2v)  o 
2G(l+v) 


(3.98f) 


£  »  0 . . . .  (3.98g) 


£  -  0 .  (3.98h) 

cpo 


e.  -  0 .  (3. 98i) 

9p 


0P9*-°U . (3-,8)> 


'  a„ . .  (3.98k) 


a80  ■  '  °u .  (3.981) 


O  *  *  0 

p4> 


(3.98m) 


<^0  -  0 .  (3.98b) 


a8p-° 


(3.98o) 
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The  solution  given  in  Eqns.  (3.98)  may  be  verified  by  an  approach  other 
than  3oussinesq's  solution  in  two  harmonic  functions.  For  a  sphere 
subjected  to  the  boundary  conditions  given  by  Eqns.  (3.94)  there  will 
be  an  infinite  number  of  possible  axes  of  symmetry  passing  through  the 
origin.  Due  to  this  symmetry  the  displacement  strain  and  stress  fields 
will  be  independent  of  the  spherical  coordinate  directions  $  and  9.  In 
view  of  this  symmetry  the  displacement  components  are  given  by 


uQ  *  u  (P) . (3.99a) 

P  p 

$  -  0 . (3.99b) 

q  -  0 . (3.99c) 


From  the  strain-displacement  relationships  given  by  Eqns.  (3.34),  the 
strain  components  are  given  by 


(3.100a) 

(3. 100b) 


’99 


(3.100c) 


’ptp 


(3. lOOd) 
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£0O  *  0 . .  (3.  lOOf ) 


From  Che  conscicucive  relationships  given  by  Sqns.  (3.37),  Che  stress 
components  are  given  by 


a  -  (X+2G)  e  +  A(e  ♦  £afl) 
PP  pp  pp  90' 


(3.101a) 


°w  '  W*2G)  *  x(£ee*  eoo) . u.ioib) 


°ee  *  a*2a)  s99  *  xtepP*  V 


(3.101c) 


O  ■  0 . (3.  lOld) 

P<p 


CT^0  -  0 .  (3.  lOle) 


aQp  -  0 .  (3.  lOlf) 


For  Che  case  when  the  stress  components  are  functions  of  the  spherical 


coordinate  p  only,  the  differencial  equations  of  equilibrium,  in  the 


absence  of  inertia  and  body  forces,  reduce  to  the  following 


dgpp  f  2ctpp  >4>  gee  m  Q 
dp  p 


(3. 102a) 


acj)(p  "  a96  *  0 


(3. 102b) 


Combining  Eqns.  (3.99),  (3.100),  (3.101)  and  (3.102)  yields  Che  follow- 
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1  up  .  2  dup  2 

,  --  r  —  — s —  —=■  a.  ■  o 

d  p  d  rt2  P 

p  P  P 


(3.103) 


The  Eqn.  (3.103)  is  a  linear,  second  order,  homogeneous  differencial 
equation  which  is  solved  for  Che  displacement  component,  u^.  Solution 
of  Eqn.  (3.103)  given  the  following  for  u^. 


u  *  A-!  P  +  — 

Pi  P 


(3.104) 


where  A^,A2  ■  constants  to  be  determined  from  the  boundary  conditions. 

Combining  Eons.  (3.100),  (3.101a)  and  (3.104),  the  stress  component, 
id  given  by 


a»0  ■  <3X  *  2  >  *1  *  -  *2 

P 


(3.105) 


The  Eqn.  (3.105)  shows  chat  in  order  for  the  stress,  ,  to  remain 

PP 

finite  at  the  origin  of .the  sphere,  Che  constant,  A2,  must  equal 
zero.  Applying  the  boundary  condition  given  by  Eqn.  (3.93a)  to  Eqn. 

(3.105),  the  constant  A^  is  given  by 


3 

A  -  u 

1  (3A+2G) 

■i 

where 

.  2VG 

A  (l-2v) 

(3.106) 


Substitution  of  Eqns.  (3.104)  and  (3.105)  into  Eqns.  (3.100)  and 
(3.101)  yields  the  following  displacement,  strain  and  stress  fields. 
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cr^g  -  0 . (3. 107n) 

^0p  -  0 . (3.  L07o) 

The  solution  given  by  2qns.  (3.107)  is  Che  seme  as  ChaC  obtained  from 
Che  Boussinesq  solution  in  Cvo  harmonic  functions.  The  Boussinesq 
solution  is  given  by  Eqns.  (3.98). 

3. 5  The  Elastic  Sphere  in  Contact  with  an  Arbitrary 
Number  of  Adjacent  Spheres 

In  this  section,  superposition  will  be  used  to  obtain  the  dis¬ 
placement,  strain  and  stress  fields  for  a  sphere  in  contact  with  an 
arbitrary  number  of  adjacent  spheres. 

In  the  remainder  of  this  Section,  tensor  notation  will  be  employed 
since  it  provides  a  convenient  way  to  express  tensor  transformations. 
The  following  rules  regarding  Censor  notation  will  be  followed. 

1.  A  superscript  denotes  a  contravariant  Censor 

2.  A  subscript  denotes  a  covariant  Censor 

3.  A  repeated  index  implies  summation  from  1  to  3. 

Superposition  of  an  arbitrary  number  of  contacts  on  a  sphere  will 
require  Che  use  of  local  and  global  coordinates  systems.  Local  coor¬ 
dinate  systems  will  be  required  for  each  contact  pair,  while  a  global 
coordinate  system  will  be  used  to  reference  the  total  solution.  Local 
coordinate  systems  set  up  for  contact  pair  tn  are  shown  in  Fig.  3.9.  a 
subscript  m  will  be  used  to  denote  the  local  coordinate  systems  for 
contact  pair  m.  The  global  coordinate  system  is  shown  in  Fig.  3.10. 


Also  shown  in  Fig.  3.10  are  the  angles  and  3n,  which  establish 
the  position  of  the  local  coordinate  system,  (x^,  ym,  z^) ,  relative 
to  the  global  coordinate  system  (x,  y,  z).  The  tensor  notation  which 
will  be  used  in  defining  the  local  and  global  coordinate  systems  is 
listed  in  Table  3.1.  Hats  (A)  and  overbars  (")  will  be  used  to  denote 
quantities  referenced  to  local  spherical  coordinate  systems  and  local 
rectangular  coordinate  systems,  respectively.  The  absence  of  these 
symbols  indicate  quantities  referenced  to  the  global  coordinate  system, 
The  tensor  components  of  the  displacement,  strain  and  stress  fields 
will  be  denoted  by  v,,  e^j  and  T^j,  respectively.  The  physical 
components  of  these  fields  will  be  denoted  by  u^,  ££j  and  T^j. 


3.5.1  Tpfle& formation  of  Displacement  Fields 

The  solution  obtained  for  the  displacement  field  resulting  from  a 
pair  of  contacts  is  referenced  to  the  local  coordinate  system,  (x-^, 

) •  To  determine  the  displacements  relative  to  the  global  coordinate 
system,  (xj_,  X£,  X3),  it  is  helpful  to  first  determine  the  displace¬ 
ments  relative  to  the  coordinate  system  (x-j_.x2.X3).  The  two  local 
coordinate  systems  are  related  through  the  mapping. 

—  A  A  A 

x^  *  x^  sinx^  cosx^ . . .  (3.108a) 

_  A.  /A 

x2  *  X1  s^nx2  3iQ*3 .  (3.108b) 

x.  *  x.  cosx_  . .  .  .  .  (3.108c) 


The  Jacobian  matrix  of  the  set  of  functions  in  Eqns.  (3.108)  is  given 
by 


3x^ 

3x2 

3i3 

3xx 

3x2 

3*3 

3x3 

«x3 

SX^ 

3^ 

3x2 

3x, 

sinx^cosx^  x^cosx^cosx^  -x^sinx^sinx^ 


sinx^inx^  x^cos^sinx^  x^sina^cosj^ 


cosx^  -x^sinx^ 


0 

(3.109) 


The  base  vectors  of  the  coordinate  system  (x^,  x^,  S3)  with  respect  to 
the  coordinate  system  (x^t  x2>  x^)  are  given  by 


o  .  -  - L  e, 

1  3*.  i 

X 


(3.110) 


where  b^  *  base  vector  in  coordinate  direction  x^. 
e^  »  unit  vector  in  coordinate  direction  x^ . 

A 

From  Eqns.  (3.108)  and  (3.109),  the  base  vectors  b^^  are  given  by 


A  A  A  —  A  A  —  A  — 

bi  *  sinx2  cosx^  e^  +  sinx^  sinx^  e2  +  cosx2  e^.  .  .  (3.111a) 


A  A  A  «  A  A  A  a  A  A  « 

b2  *  x^  cosx2  cosx^  e^  +  x^  cosx2  sinx^  e2  -  x^  sinx2  e^ 


b-  •  -x.  sinx_  sinx_  e.  +  x, sinx-  cosx. 


(3.111b) 

(3.111c) 
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The  metric  tensor  is  determined  from  the  dot  products  of  the  base 
vector  a3  follows 


«y  -  si  '  SJ 


(3.112) 


where  g^_.  ■  component  of  the  metric  tensor 


From  Eqn.  (3.112),  the  metric  tensor  is  given  by 


0  (xi)2 


0  (xj_)2  sin^  X2 


03.113) 


The  tensorial  components  of  the  displacement  field  in  the  coordinate 
system  (x^,  ,  x^)  are  given  by 


1  /  A  A  /A.  \ 

u  (j^.x^x  ) 


v1(x1,x2,x3) 


,  no  sum  on  i . (3.114) 


where 


v1  *  components  of  the  displacement  tensor 

A  . 

u1  *  physical  components  of  the  displacement 


tensor 


To  determine  the  displacement  components  vlt(x]_1X2,X3) ,  the 
transformation  law  for  a  contravar iant  tensor  of  order  1  is  used.  This 
transformation  law  is 


V  (x1,x2,x3)  =  V  (x1,x2,x3)  — 7“ . (3.115) 

3x_. 

From  Eqns.  (3.109),  (3.113)  and  (3.115),  Che  displacement  components 
-k  -  -  -  . 

v  (x2,X2,x3)  are  given  by 


(v>  »  [Tl]T  {u} 


(3.116) 


The  vectors  and  matrix  appearing  in  Eqn.  (3.116)  are  given  by 


{v} 


vl 


-2 

v 


vJ 


/  \ 

u. 


\  u3  i 


(3.117a) 


{u} 


t  ~  > 

“1 


< 


u„ 


\  .  (3.117b) 


'  u3  ' 


sinx^ 

cosx^ 

sinx^ 

sinx  2 

cosx. 

IT  I  - 

cosx2 

cosx^ 

cosx^ 

sinx^ 

-sinx. 

-sinx. 


cosx. 


0  J 


(3.117c) 


In  Eqns.  (3.117)  the  physical  components  u^  are  used  in  place  of  the 
tensor  components  v\  In  a  rectangular  coordinate  system,  no  distinction 
is  made  between  the  contravariant,  covariant  and  physical  components  of 


a  tensor. 
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The  global  coordinates,  (x^,  x^,  x^)  are  related  to  the  local  coor¬ 


dinates,  (x^,  x^  x^) ,  by  the  mapping 


x.  *  cos 3  cos  jL  +  -sin\|;  -  sinB  cosip  x-  .  (3.118a) 

1  m  ml  m  m  m  j 

X,  ■  cos3m  sin  \jj  x.  +  cosii  x„  +  sinB  sini|>  x.  .  (3.118b) 

*  m  mi  Tni  2  m  m  J 

X  -  -sin6_  X  -  cos3„  X . (3.118c) 

J  ml  m  3 

The  Jacobian  matrix  [J]  of  the  functions  in  Eqns.  (3.118)  is  given  by 


cos3 

m 

cosij; 

m 

-sinij/ 

m 

sinB 

m 

COSli) 

cos  6 

m 

sinijj 

m 

costy 

m 

sinB 

m 

sinijj 

-sinB 

m 

0 

cosB 

m 

The  displacement  field  in  the  global  coordinate  system  may  be  determined 
from  the  transformation  law  applicable  to  a  contravariant  tensor  of  order 
one.  This  transformation  is  given  by 

.  .  3x. 

v  (X;L,x2,s3)  -  ^>(x1,x2,x3)  . (3.120) 

j 

From  Eqns.  (3.119),  the  displacement  components  v^x^.x^x^  are  given 
by 


(v>  -  [T  ]  T  (v) 


(3.121) 


I 


The  vectors  and  matrix  appearing  in  Eqn.  (3.121)  are  given  below 


(3.122a) 


(3. 1.2b) 


kl 


COSg  COS  til 

m  Tm 

-sin’l' 

m 

sin  3  cos<P 
01 

cos  3  sin'ii 
m  m 

COSlj/ 

m 

sinS  sin^i 
01 

-sinS 

m 

0 

cos  0 

01 

In  Eqns.  (3.122),  the  physical  components  of  displacement  are  used 
since  no  distinction  is  made  between  the  tensor  and  physical  components 
in  a  rectangular  coordinate  system.  When  the  displacement  components 
u^(x^,X2»x3^  are  known,  the  global  displacements  u^(x^,X2,X2^  may 
be  determined  from  Eqns.  (3.116)  and  3.121).  The  displacement  compo¬ 
nents  u^(x^,X2,X2)  are  given  by 


J  u[  -  (Tg]T  [TlJT  |u [ 


(3.123) 


3.5.2  Transformation  of  Strain  Fields 


The  solution  to  the  strain  field  resulting  from  three  spheres  in 


contact  along  an  axis  of  symmetry  is  referenced  to  the  coordinate 


system  (x^.Xj.j^).  To  determine  the  strain  field  in  the  global  coor¬ 
dinate  system,  (x  ,x  ,x,),  it  is  helpful  to  first  determine  the  strain 
field  in  the  local  coordinate  system,  (x^^  .x^) .  The  relation  between 


the  coordinate  systems  (x^.x^.x^)  and  (x^.x^.x^)  is  given  by  Eqns. 


(3.108).  The  Jacobian  matrix  [j]  of  the  functions  defined  in  Eqns. 


(3.108)  is  given  by  Eqns.  (3.109).  The  metric  tensor  of  the  coordinate 


system  (x^.x^.x^)  with  respect  to  the  coordinate  system 


is  given  by  Eqn.  (3.113).  The  tensorial  components  of  the  strain  field 


in  the  coordinate  system  (x^.x^.x^)  are  given  by 


a1]  f  ^  /s,  vs  x 

e  J(x1>x2>x3) 


,  no  som  on  i  or  j . (3.124) 


8ii  y  8jj 


where 


components  of  the  strain  tensor 


*  physical  components  of  the  strain  tensor. 


To  determine  the  strain  field  e^  (xl»x2,x3^ ,  the  transformation  law  for 


a  contravariant  tensor  of  order  two  is  used.  This  transformation  law 


-ij,-  -  -  .  Jdi,.  -  -  x  ifi  ifl 

e  (X3_.X2.X3)  *  e  vX3_.X2.X3)  ^ 


(3.125) 


From  Eqns.  (3.109),  (3.124)  and  (3.125),  the  components  of  the  strain 


tensor,  e^(x,,x„x J,  are  given  by 
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i\r  le^i  [tl] 


(3.126) 


The  transformation  matrix  [TjJ  is  given  by  Eqn.  (3.117c).  The 

ij 


matrices  [e  ]  and  [£  ]  are  given  by 
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(3.127a) 
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£l  3 
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2  1 
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(3.127b) 


In  Eqn.  (3.127a)  the  physical  components  may  be  used  in  place  of  the 
tensor  components,  since  the  strain  tensor,  [e^],  is  referenced  to  the 
rectangular  coordinate  system,  (x^.Xj.x^). 

The  relationships  between  the  global  coordinates,  (x^,X2,X2)  and 
the  local  coordinates,  (x^.Xj.x^),  are  given  by  Eqns.  (3.118).  The 
Jacobian  matrix  [J]  of  the  functions  given  by  Eqns.  (3.118)  is  given  by 
Eqn.  (3.119).  To  determine  the  strain  field  in  the  global  coordinate 
system,  the  following  transformation  law  for  a  second  order  tensor  will 
be  used. 


—  il  _  ]r0  —  —  — 

e  ■'(x^^.x^)  »  e*cx'(x^,X2>Xj )  — — =■  — ~J- 

5x,  9x 
k  e 


(3.128) 


From  Eqns.  (3.119)  and  (3.128),  the  strain  field  e  J (x1 ,X2 .x^)  is  given 
by 
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teij]  »  [TG]T  [iij]  [TrJ 


(3.129) 


The  matrices  [Tq]  and  [e  are  given  by  Eqn.  (3.122c)  and  Eqn. 
(3.127a),  respectively.  The  strain  field  [e1J]  is  given  by 
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31  32  33 
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.  (3.130) 


In  Eqn.  (3.130),  the  physical  components  of  the  strain  tensor  may  be 
used  as  the  tensor  components  since  the  strain  field  is  referenced  to 
the  rectangular  coordinate  system  (x^.x^.x^).  Combining  Eqns.  (3.126) 
and  (3.129),  the  strain  field  [e^l  is  given  by 


[e1J]  -  [TG1T  [Tl]T  [«  J  [Tl]  [Ic] 


(3.131) 


The  matrices  [TLJ  ,  [Tq]  and  [e^]  are  given  by  Eqns.  (3.117c), 
(3.122c)  and  3.127b),  respectively. 


3.5.3  Transformation  of  Stress  Fields 

The  stress  fields  in  the  rectangular  coordinate  systems  (x^.x^.x^) 
and  (x^^^.a^)  are  determined  in  the  same  manner  as  the  strain  fields. 
Both  the  stress  and  strain  tensors  are  second  order  and  vill  transform 
in  the  same  manner.  Therefore,  the  results  determined  for  the  strain 
fields  in  section  3.5.2  may  be  used  to  determine  the  stress  fields  in 
the  rectangular  coordinate  systems  (x^.x^.x^)  andte^.x^.x^).  From 


these  results,  the  stress  fields  are  given  by 


[Tlj]  -  [Tl]T  [o  ]  [Tl1 . (3.132) 

[tij]  -  [TG]T  [?ij]  [Tg]  . (3.133) 


The  matrices  (T^I  and  [TqJ  are  given  by  Eqn.  (3.117c)  and  Eqn. 
(3.122c),  respectively.  The  matrices  [x1J]  and  [x1^]  are  given  by 
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(3.134) 


(3.135) 


In  Eqns.  (3.134)  and  (3.135),  the  physical  components  of  the  stress 
tensor  may  be  used  for  the  tensor  components,  since  the  stress  fields 
are  referenced  to  rectangular  coordinate  systems.  Combining  Eqn. 
(3.132)  and  Eqn.  (3.133),  the  stress  field  in  the  global  coordinates 
(xj_.x2.X3)  is  given  by 


[Tij]  -  [tg]t  [Tl]T  [5tj]  [tl]  [tg] 


(3.136) 


The  matrices  [T^]  and  [T^]  are  given  by  Eqn.  (3.117c)  and  Eqn.  (3.122c), 
respectively. 
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3.5.4  Solution  for  a  Sohere  with  an  Arbitrary  Nuusber  of  Contacts 


The  results  of  Sections  3.3.1,  3.5.2  and  3.5.3  will  be  used  to 
obtain  the  displacement,  strain  and  stress  fields  for  a  sphere  subject 
to  an  arbitrary  number  of  contacts  with  adjoining  spheres.  Restric¬ 
tions  which  are  imposed  on  these  solutions  are  listed  below. 

1.  The  contacts  on  the  sphere  appear  in  pairs  along  an  axis  of 
3ymmetry  as  shown  in  Fig.  3. 9 

2.  The  forces  transmitted  through  the  contacts  are  such  that  the 
sphere  is  in  static  equilibrium. 

Subsections  3.51,  3.52  and  3.53  provide  the  displacement,  strain  and 
stress  fields  referenced  to  a  global  coordinate  system.  These  fields 
were  developed  for  the  case  of  one  pair  of  contacts  along  an  axis  of 
symmetry  as  appearing  in  Fig.  3.9  If  a  sphere  has  a  number  of  these 
contact  pairs,  the  displacement,  strain  and  stress  fields  resulting 
from  these  contacts  may  all  be  referenced  to  the  global  coordinate 
system.  The  displacement,  strain  and  stress  fields  resulting  from  all 
contacts  are  then  found  by  adding  the  individual  fields.  Superposition 
of  these  individual  fields  yields  the  following  total  displacement, 
total  strain  and  total  stress  fields. 
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(3. 137a) 


(3. 137b) 
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where  nc  *  total  number  of  contacts  on  the  sphere. 

In  Eqns.  (3.137),  the  subscript  T  is  used  to  denote  the  total  displace¬ 
ment,  strain  and  stress  fields.  The  subscript  m  denotes  quantities 
resulting  from  the  pair  of  contacts,  m. 

The  quantities  appearing  in  Eqns.  (3.137)  may  be  put  in  terms  of 
th*  local  coordinate  system  and  global  coordinate  system  shown  in  Fig. 
3.9  and  Fig.  3.10  respectively.  The  total  displacement,  strain  and 
stress  fields  may  be  written  as  follows 
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(3.138b) 


(3.138c) 


In  Eqns.  (3.138),  the  global  coordinate  system  (x,  y,  z)  is  referenced. 
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Che  displacement,  strain  and  stress  fields  resulting  from  contact  m  are 
given  by 
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(3.139a) 


(3.139b) 


(3.139c) 


The  fields  lu  ]B  and  ]B  are  determined  from  the 

results  given  in  Section  3.4.  These  fields  are  referenced  to  the 
spherical  coordinate  system  (pB,  0B>  0B)  shown  in  Fig.  3.9  .  The 
transformation  matrix,  [ Tq ] ,  is  given  by 


cos  B_  COS<j; 

m  m 
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sinS  cos0 
m  m 

cos  3  cos0 
m  m 

COS0 
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sing  sinui 
a  m 

-sin 8^ 
m 

0 

cos  g 
m 

(3.140) 


The  angles  8^  and  0  define  the  position  of  the  local  rectangular  coor¬ 
dinate  system,  ( xm . ym . )  relative  to  the  global  coordinate  system, 


m  ’  «  *  »  -  *  -  ■  •*  ■  > 


(.x,  y,  z).  These  angles  are  shown  in  Fig.  3.10.  The  trans  format  ion 
matrix,  [TjJ  ,  is  given  by 


sin9  cos  @ 
m  m 

cos<J>  COS0 
m  m 

-sin@ 


sind>  sin3 
m  m 

COS0  cos9 
m  m 


-sind> 


(3.141) 


The  angles  9n  and  0a  are  shown  in  Fig.  3.  9  .  These  angles  are  part 
of  the  spherical  coordinate  system,  (f^,,  4^,  9m).  They  define  a 
point  in  the  sphere  relative  to  the  local  coordiante  system,  (x^  ym,  z 
for  contact  m. 

In  order  to  evaluate  the  total  displacement,  strain  and  stress 

fields  at  any  point  (x,  y,  z),  the  local  coordinates  must  be  determined 

from  the  global  coordinates.  This  is  necessary  since  the  solutions 

provided  in  Section  3.4  are  in  terms  of  the  local  coordinates  (p_,  9  , 

“  -m 

The  following  equations  are  used  to  determine  the  local  coordinates 
(*m»  7m*  *m^  and  (Pm>  $m»  em^  when  tbe  global  coordinates  (x,  y,  z)  are 
known. 


cos  S  cosip  cosB  sinti;  -sing 
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,  m  m 


COS9 
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sinS  sin9  cos  6 
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V2  2  2 

x  +  y  +  z . (3.143a) 


L  2  .  2 

x  +  y 

m  ym 


(3.143b) 
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(3.143c) 


(3.143d) 


(3.143e) 


From  Eqns.  (3.142)  and  (3.143)  the  total  displacement,  strain  and  stress 
fields  may  be  evaluated  for  any  point  (x,  y,  z) . 
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CHAPTER  4 

EFFECTIVE  MODULI  OF  AN  IDEALIZED  SOLID  -  FLUID  SYSTEM 


4.1  Effective  Moauli 


.The  effective  moduli  of  heterogeneous  materials  are  those  moduli 
reflecting  the  average  stress-strain  properties  of  the  materials.  These 
moduli  taKe  into  account  the  properties  and  the  geometry  of  all  the  phases  and 
tneir  interaction.  The  effective  moduli  are  determined  by  considering  a  small 
but  representative  sample  of  the  heterogeneous  materials.  The  size  of  the 
representative  sample  is  chosen  so  that  the  behavior  of  this  sample  does  not 
cnange  when  the  sample  size  is  increased.  Therefore,  the  effective  moduli 
determined  for  the  sample  may  be  used  to  represent  the  entire  heterogeneous 
material.  The  two  approaches  that  are  used  to  determine  the  effective  moduli 
are  volumetric  averaging  and  energy  methods. 

4.1.1  Volumetric  Averaging  Approach 

The  volumetric  averaging  approach  derives  effective  material  properties 
by  considering  a  representative  volume  element  of  the  material.  The  stress  or 
displacement  fields  in  this  representative  volume  element  are  macroscopically 
homogeneous.  The  volume  averaged  stress  is  defined  as 


<°,  -  v  f  o,  ,  d V 


ij 


ij 


(4.1  ) 


wnere  <&  >  “  volume  averaged  stress  field 


a  -  stress  field  occurring  in  the  different  phases  of  the  material 

•  *J 


V  -  total  volume  of  representative  sample 
In  Eqn.  (4.1),  Cartesian  tensor  notation  has  been  used.  This  notation  will  be 
used  periodically,  and  the  rule  that  a  repeated  tensor  implies  summation  is 
applicable.  The  volume  averaged  strain  is  defined  as 
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<£U>  =  V  i  eij  aV 


(4.2) 


<£.  >  =*  volume  averaged  strain  field 


e.  .  -  strain  field  occurring  in  the  different  phases  of  tne 
^  d 


materials 


Both  of  tne  integrals  appearing  in  Eqns.  (4.1)  and  (4.2)  are  taken  over  tne 


representative  volume  of  the  material.  The  effective  moduli  are  defined 


througn  tne  following  equation. 


<a  .>  =*  C  ,  ,  <e.  ,> 
ij  ljki  kJt 


(4.3) 


wnere 


effective  moduli  tensor 


The  Eqn.  (4.3)  is  the  general  anisotropic  form  of  the  linear,  volume  averaged, 


stress-strain  relations.  In  order  to  use  Eqn.  (4.3)  to  determine  the 


effective  moduli  of  a  heterogeneous  material,  the  stress  and  strain  fields 


occurring  in  all  phases  must  be  known  and  the  integrals  in  Eqns.  (4.1)  and 


(4.2)  performed. 


4.1.2.  Energy  Methods 


When  using  an  energy  approach,  the  effective  moduli  are  determined 


througn  energy  equivalence,  a  representative  sample  of  the  heterogeneous 


material  is  considered  to  be  subjected  to  a  macroscopicaliy  homogeneous  stress 


or  deformation  field.  The  sum  of  the  strain  energies  for  all  the  phases  of 


the  neterogeneous  material  is  then  made  equivalent  to  that  occurring  in  a 


Homogeneous  material  with  moduli,  This  relationship  is  given  by 


eij‘Vv '  <e1J>  7 


(4.4) 


•M 


mumm 


wnere 


stress  rieid 


=  strain  fieia 

V  =  total  volume  of  representative  sample. 

The  volume  averaged  stress,  <0^^,  and  the  averaged  volume  strain,  <e1j>.  are 
determined  from  Eqns.  (4.1)  and  (4.2),  respectively.  Substituting  Eqn.  (4.3) 
into  Eqn.  (4.4)  yields  tne  following 


J  0  .  £  .  dV 

v 


C  <£  >  <£,  ,  > 

ljkfl,  ij  k  i 


(4.5) 


Tne  integral  appearing  in  Eqn.  (4.5)  is  taken  over  the  entire  representative 
volume.  As  in  the  volumetric  averaging  approach,  the  stress  fields,  o. .,  and 
the  strain  fields,  e. .,  must  be  known  for  the  different  phases  of  the 

*  *J 

neterogeneous  material  before  Eqn.  (4.5)  may  be  evaluated. 

4.2  Volumetric  Averaging  Approach  to  Determining  the  Effective  Bulk  Modulus 

of  an  Idealized  Solid-Fluid  System 

The  solid-fluid  system  to  be  considered  is  a  heterogeneous  material 
consisting  of  a  number  of  solid  spheres  in  contact,  surrounded  by  a  fluid 
pnase.  Sucn  a  system  is  snown  in  rig.  4.1.  The  following  assumptions  are 
used  in  the  determination  of  the  effective  bulk  modulus. 

1 .  The  spneres  contained  in  the  system  consist  of  isotropic,  linearly 
elastic  materials. 

2.  The  spheres  contained  in  tne  system  are  m  static  equilibrium. 

3.  The  contacts  on  a  particular  sphere  occur  in  pairs,  directed  along  an 
axis  passing  tnrougn  the  spnere. 

4.  The  surface  displacements  or  surface  tractions  resulting  from  a  pair 
of  contacts  are  axisymmetric  witn  respect  to  an  axis  passing  tnrougn 
tne  center  of  tne  contact  region. 


5.  The  solid-fluid  system  benaves  as  an  isotropic,  linearly  elastic 
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material  macrascopica.uiy . 

o.  The  voia  space  between  spneres  in  the  system  is  completely  filled 
with  only  one  type  of  fluid. 

For  an  isotropic,  lineany  elastic  material,  only  two  material  constants  are 
needed  to  relate  tne  force  field  to  tne  displacement  field.  For  this  case, 
Sqn.  (4.3)  reduces  to 


<o  >  -  2G  <e  >  +  A  <e  >  6 
ij  ij  kk  ij 


(4.6) 


wnere  <aij>  ”  volume  averaged  stress  tensor 
<e  >  *  volume  averaged  strain  tensor 
d,A  »  effective  Lame  constants 


6  .  -  Kronecker  delta  -  i°  1  *  J 
ij  i  ,  i  -  J 


Anotner  relationship  may  be  obtained  by  contracting  the  i  and  j  indices  in 
Eqn.  (4.6).  This  contracted  form  is  given  by 


<akk>  *  3  K  <ekk> 


(4.7) 


where 


3  2  -  n 

K  -  J  G  +  A 


tne  term  K,  appearing  in  Eqn.  (4.7),  is  the  effective  bulk  modulus.  The 
quantities  <a^>  and  <£ki<>>  appearing  in  Eqn.  (4.7),  are  determined  from  Eqn. 
(4.1)  and  (4.2) ,  respectively.  These  quantities  are  given  by 


<akk>  V  ^  akk 


(4.8a) 


<e  >  »  —  f  p  aV 
kK  ^  V  kk 


(4.8b) 
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a  =  trace  of  tne  stress  tensor 

KK 


trace  of  tne  strain  tensor 


V  -  representative  volume 


Trie  integrals  appearing  in  Eqns.  (4.3)  are  taKen  over  tne  volumes  of  Doth  the 
solid  and  the  fluid  pnases  of  the  heterogeneous  material.  Due  to  this,  the 
integrals  appearing  in  Eqns.  (4.8)  may  be  rewritten  as  follows. 


<a.  ,  > 
kk 


<e,  ,  > 
kk 


1  s  f 

V  L  /  °kk  d V  -  ;  akK  aV  j 


V  L  i  ekk  dV  +  fv  £kk  dV  j 

3  f 


(4.9a) 


(4.9b) 


where  Vg  =  volume  of  tne  solid  pnase  of  tne  neterogeneous  material 
-  volume  of  the  liquid  phase  of  the  heterogeneous  material 
The  superscripts  s  and  f  appearing  in  Eqns.  (4.9)  are  used  to  denote  tne  solid 
and  fluid  phases,  respectively.  Combining  Eqns.  (4.7)  and  (4.9),  the 
following  expressions  for  the  effective  bulk  modulus,  K,  is  obtained 


?  L  /  okk  dV  «•  /  aRk  dV  j 

_ _s _ _f _ 

V  !•  *  „  ekk  dV  +  Sck  dV  ■ 


(4.10) 


To  evaluate  Eqn.  (4.10;  a  representati ve  volume  of  tne  solid-fluid  system  is 
considered.  This  representative  volume  is  taken  as  tne  smallest  possible 
sample  wrucn  still  exmbits  the  same  benavior  as  the  entire  system.  Due  to 
this,  tne  effective  bulK  modulus,  tf,  found  by  performing  the  integrals  m  E^n. 
(4.10)  over  this  representative  volume,  will  be  applicable  to  the  entire 
system.  In  determining  tne  effective  bulk  modulus,  K,  the  stress  and  the 
displacement  fields  occurring  in  the  representative  volume  will  be  considered 


V 


£ 


I 


8 


:S 


£ 


& 


a 


Ni 


>7 

f.;1 


i  v, 


co  be  macroscopically  nomogeneous. 

Tne  representative  volume  wiix  contain  N  soiia  spneres.  The  apneres  are 
considered  to  be  isotropic,  linearly  elastic  materials.  They  may  nave 
different  raan  and  material  constants.  These  spneres  constitute  tne  solid 
phase  of  the  heterogeneous  material  shown  in  Fig.  4.1.  In  Eqn.  (4.10).  the 
integrals  taken  over  the  soiid  pnase  of  the  representative  volume  are -written 
as  follows. 


where 


i  si;  s 

V  °kk  aV  *  V  4  *  A  Kk)x  aV 

s  s 

1  N 

V  l  ekk  dV  "  V  i.  ;ui  (ekk}l  dV 

V  1*1  V 

S  3 

s 

(akk}i  “  trace  of  stress  tensor  for  sphere  i. 

(£Wi  *  Crace  of  cl»e  strain  tensor  for  sphere  i. 


(4.11a) 


(4.11b) 


In  Eqns.  (4.11)  the  integrals  appearing  in  the  summations  are  taken  over  the 
volume  of  spnere  i,  V  .  Due  to  the  first  four  of  tne  assumptions  made 
concerning  the  solid-fluid  system,  the  results  in  Section  3.3  may  be  used  to 
determine  (okl<)  and  ( e^)  .  From  Section  3.3,  the  solution.  [S],  to  the 
displacement,  strain  or  stress  fields  for  a  sphere  undergoing  axisymmetnc 
surface  tractions  or  axisymmetnc  surface  displacements  are  given  by 
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£ 
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CO  00 

[3]  ■  \  Sn-rA-„-,J  *  ; 

n-1  n-o 

wnere  [S]  «  solution  to  tne  displacement,  strain,  or  stress  field 

L A _n_ ^  J  ■  component  solution  to  the  displacement,  strain  or  stress 
f  leld 

C ^_n_2 J  "  component  solution  to  tne  displacement,  strain  or  stress 


field 
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=  constant  of  superposition 
*  constant  of  superposition 

Tne  component  solutions,  [A_n-1 j  ana  [ 3_n_^ J  >  can  be  determined  from  Eqns. 
(3.46)  ana  (3.53).  respectively.  The  constants  of  superpositions,  a_n_-,  ana 
d_n_2»  are  aetermineC  from  tne  specifiea  Doundary  conaitions  on  the  sphere. 
Evaluating  (a  )  and  (e  )  corresponding  to  the  component  solution  [A_  . ] 

iK K  X  KK  X  n  ■  * 

gives  the  following 


-in-1 


“W1 


(el> 


-nnl 


0 


(4.13a) 


(4.13b) 


In  Eqns.  (4.13),  the  superscript  A_n-1  is  used  to  denote  quantities  which  are 

3  S 

determinea  from  tne  component  solution  [A  _1].  Evaluating  (a^)^  and  C  G-kt<) 
corresponding  to  tne  component  solution  [B_n_2]  yields  the  following 


(o  *)  -  2Cn+l  )(2n«-3)0+v  )  pnPn  (4.14a) 

B-n-2  (1-2v. ) 

(ekk)i  »  (n+1 )  (2n+3)  1  pnpn  (4.14b) 

wnere  v  -  Poisson's  ratio  for  sphere -i. 

G  -  shear  modulus  of  sphere  l. 
p  -  distance  from  tne  origin  to  a  point  in  sphere  l. 

In  Eqns.  (4.14),  tne  superscript  3_n_2  is  used  to  denote  quantities  determined 
from  tne  component  solution  Combining  Eqns.  (4.12)  through  (4.14) 

3  3 

yields  the  following  for  (o  )  and  (e  )  ,  when  there  is  one  pair  of  contacts 

KK  X  K  K  x 

on  sphere  i. 


(c^).  -  2  S  (b.  ,).(n+l)(2n+3)(l-vJo“  P. 


-n-2  i 


( ef . )  Z  (b  _).(n+l)(2n+3)(l-2v.)  pD  P  . 

kk  G.  ,  -n-2  i  i  n 

i  n-1 


(4.15a) 


(4.15b) 


The  subscript  i  is  used  on  the  constant  of  superposition,  b  since 

its  value  will  depend  on  sphere  i  and  the  spheres  making  contact  with 
sphere  i.  This  dependence  is  shown  in  Section  3.1  where  the  Hertz 
contact  problem  is  discussed.  The  quantities  (a,?),  and  (  e®  )  .  are 

l&iC  X  1C1C  X 

invariant  with  respect  to  coordinate  directions.  Therefore,  for  a 

s  s 

sphere  with  M  contacts,  ( O^^)  and  (  e^)  are  given  by 


(O  *)  . 

kk  l 


n  r  .  2  ki 

.  *  2  I  (b  ,).+  2  (b  .  I  (n+l)(2n+3)(l+v- )  P  Pn 

1  n-lL  "“‘2  1  k-1  -°-2  1J  1  n 


(4.16a) 


,  S  X  1 

(e. , )  .  ■ 

kk  i  G. 
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[<5-n-2>i+ 


(n+l)(2n+3)(l-2VjL)  D  Pn 


(4.16b) 


where  Mj_  *  number  of  contacts  on  sphere  i. 


In  Eqns.  (4.16),  the  superscript  k  has  been  used  on  the  constant  of 
superposition,  b_^_2 »  since  this  quantity  will  change  for  the  different 
contacts,  occurring  on  a  single  sphere.  The  constant  of  superposition, 
b_Q_2 i  results  from  the  application  of  a  radial  pressure  to  the  surface 
of  sphere  i.  Substituting  Eqns.  (4.16)  into  Eqns.  (4.11)  and  perform¬ 
ing  the  indicated  integration  yields  the  following 


V  akk  dV 
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(4. 17a) 


v  4 dv 
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S  (b-2>i 
k-1 


(l-2v.)R7 

l  l 


(4.17b) 


where  -  Che  radius  of  sphere  i. 


To  determine  the  integrals  involving  cr^  and  e^,  which  appear  in  Eqn. 


(4.10),  the  fluid  pressure  in  the  representative  volume  is  denoted  as 

•  f  .  Pf 
kk  and  kk 


a  .  The  quantities  O'.  ^  and  are  related  to  ou  by 


a. .  »  3  a 

kk  u 


(4.18a) 


£  -  -  r 

'kk  I 

*  n 


au 


dP 


Kf(P,T)  f 


(4.18b) 


where  P  -  fluid  pressure 

Kj  -  bulk  modulus  of  the  flued 
T  -  temperature 

In  Eqns.  (a..  18)  the  fluid  pressure,  Ou,  is  interpreted  to  be  a 


compressive  stress  and  therefore  its  magnitude  is  negative.  This 
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pressure  is  also  a  gauge  pressure  meaning  chat  if  the  fluid  pressure  is 

atmospheric,  cu  is  equal  to  zero.  The  Eqn.  (4.18)  states  that  in 

general,  the  bulk  modulus  of  a  fluid,  Kf,  is  not  constant.  It  may 

vary  with  the  pressure  and  temperature.  The  lower  limit  of  the 

integral  appearing  in  Eqn.  (4.18b)  indicates  that  the  fluid  will  be 

considered  to  be  unstrained  at  zero  gauge  pressure.  Henceforth,  the 

f 

volumetric  strain  of  the  fluid,  C^k*  be  denoted  by  ef.  Substi¬ 

tution  of  Eqns.  (4.18)  into  the  integrals  given  in  Eqn.  (4.10)  yields 
the  following. 


where 


(4.19a) 


(4.19b) 


V  *  representative  volume 

Vf  *  volume  of  fluid  present  in  the  representative  volume 

ou  "  fluid  pressure  in  representative  volume 
ef  *  volumetric  strain  of  the  fluid  in  the  representative 
vo lume 


The  term,  Cf,  appearing  in  Eqns.  (4.19)  is  the  volume  fraction  of  the 
voids  containing  fluid,  in  the  solid-fluid  system.  Combining  Eqns. 


(4.10),  (4.17)  and  (4.18)  given  Che  following  expression  for  the  effec¬ 
tive  builk  modulus  of  the  solid-fluid  system  under  consideration. 


141 


«• 

t 

.i 

i 

1 

‘i 

1 

»< 

'i 

! 

I. 

t 


1  ‘  4il 

f  L  Gi 


(b  +  I  (b  (1+v.)  R3  +  3  C-  5 

-l  i  -2  ij  ii  f  u 


(b_2).  +  Z  (b_2)^  (l-2v.)  R3  +  Cf  ef 


(4.20) 


The  Eqn.  (4.20)  gives  the  effective  bulk  modulus,  K,  of  a  solid-fluid 


system  composed  of  a  number  of  isotropic,  linearly  elastic  spheres  in 
contact  and  surrounded  by  a  fluid  phase.  The  assumptions  made  in 
arriving  at  Eqn.  (4.20)  are  listed  at  the  beginning  of  this  section. 


Most  of  these  assumptions  were  made  so  that  the  results  from  Section 
3.2  could  be  used.  To  evaluate  Eqn.  (4.20),  a  knowledge  of  the  sizes, 
material  properties,  number  of  contacts,  location  of  contacts,  forces 
on  contacts  and  the  materials  in  contact  have  to  be  known.  It  is 


unlikely  that  these  variables  would  be  known  in  a  deterministic  form. 
Therefore,  the  evaluation  of  Eqn.  (3.20)  would  require  statistical 


: 


data.  In  the  next  subsection  some  simplications  will  be  made 
concerning  the  solid-fluid  system,  which  will  allow  an  evaluation  of 
Eqn.  (4.20). 

?rom  Eqn.  (4.20)  some  observations  can  be  made  concerning  the 
effective  bulk  modulus,  K.  For  the  case  when  the  value  of  ou  is  zerc 
gauge  pressure,  Eqn.  (4.20)  shows  that  the  effective  bulk  modulus 
depends  only  on  the  material  properties  and  geometries  of  the  spheres 
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included  in  the  solid-fluid  system.  The  same  results  is  obtained  by 
having  the  volume  fraction  of  the  voids  in  the  system,  Cf,  equal  to 
zero.  Although  the  two  cases  listed  above  give  the  same  result,  che 
manner  in  which  these  results  were  obtained  are  different.  For  the 
case  of  zero  gauge  pressure  the  solution  for  a  sphere  in  contact  with 
two  adjacent  spheres  must  be  considered.  This  solution  allows  Che 
sphere  in  contact  to  deform  freely  into  the  surrounding  void  space  for 
conditions  of  zero  gauge  pressure.  The  effect  of  this  is  to  consider 
chat  the  air  in  the  voids  does  not  contribute  to  the  stiffness  of  the 
solution.  The  Eqn.  (4.20)  would  yield  a  different  result  if  was 
an  absolute  pressure.  Under  these  conditions  the  term  ef  would  still 
be  zero  and  che  contribution  of  ou  in  che  numerator  of  Eqn.  (4.20) 
would  be  negligible.  For  this  case  when  the  volume  of  voids  is  equal 
to  zero,  the  system  would  only  contain  the  solid  phase  so  that  the 
effective  bulk  modulus,  K,  would  depend  only  on  the  material  properties 
and  geometries  of  spheres  in  the  system.  In  this  case  the  system  would 
consist  of  only  spheres  such  that  all  the  void  space  being  occupied  by 
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a  sphere  of  proper  size.  Another  case  to  consider  is  when  |  ou  |  >  0, 
and  the  fluid  phase  is  incompressible.  For  this  case  the  volumetric 

strain  of  the  fluid  phase,  eg,  is  zero.  To  study  the  effect  of  the 
incompressible  fluid  phase,  it  is  recognized  that  the  terms  appearing 

in  the  numerator  and  the  denominator  of  Eqn.  (4.20),  will  have  negative 

values  when  the  solid-fluid  system  is  in  compression.  Due  to  this, 

Eqn.  (4.20)  shown  that  Cl. »  effective  bulk  modulus,  K,  increases  as  j  J  i 

u' 

increases,  for  a  constant  volume  of  solids  and  voids.  This  shows  that 
a  non-zero  pressure  in  an  incompressible  fluid  phase  serves  as  a 
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restraint  against  volume  de f ormat ion.  This  is  reasonable  since  an 
incompress ib le  fluid  phase  should  increase  the  stiffness  of  the  system 
and  the  fluid  pressure,  J  u ,  should  exist  for  such  a  fluid.  For  the 
limiting  case  when  the  solid  phase  is  absent  in  the  system,  Eqn.  (4.20) 
shows  that  the  effective  bulk  modulus,  K,  becomes  infinite.  This  is 
reasonable  since  the  system  contains  only  an  incompressible  fluid 
phase.  For  a  compressible  fluid  phase  in  the  solid-fluid  system,  the 
volumetric  strain  of  Che  fluid  phase,  ef,  will  be  non-zero.  In  the 
limiting  case  when  the  system  contains  only  a  compressible  fluid,  the 
effective  bulk  modulus,  K,  in  Eqn.  (4.20)  is  a  function  of  the  bulk 
modulus  of  the  fluid.  This  dependence  is  apparent  from  Eqn.  (4.18b). 
For  the  case  of  a  fluid  with  a  constant  bulk  modulus,  Eqns.  (4.18b)  and 
(4.20)  show  that  the  effective  bulk  modulus,  K,  is  equal  to  the  bulk 
modulus  of  the  fluid. 

The  effective  bulk  modulus  given  by  Eqn.  (4.20)  was  determined  by 
considering  a  small  representative  sample  of  the  solid-fluid  system. 

The  macroscopic  system  will  be  made  up  of  a  number  of  these  small 
representative  samples.  If  the  displacement  or  the  stress  fields  vary 
with  Che  position  in  the  macroscopic  sample,  the  effective  bulk  modulus 

will  vary  with  position.  This  occurs  due  to  the  fact  that  the  constant 
of  superposition,  b_2,  is  dependent  on  the  boundary  conditions 

present  on  the  contacts  between  spheres.  These  boundary  conditions 
change  with  the  force  or  displacement,  on  the  contact. 

4.2.1  Effective  Bulk  Modulus  of  a  Solid-Fluid  System  Consisting 
of  E^ual  Spheres  Arranged  in  Ideal  Packing  Conf igurat ions , 
and  Surrounded  by  a  Fluid. 

Restriction;-  concerning  the  geometry  and  material  properties  of 
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Che  solid  phase  present  in  Che  soLid-fluid  system  will  be  made.  These 
restrictions  will  allow  Eqn.  (4.20)  to  be  evaluated  in  the  absence  of 
statistical  data.  The  restrictions  on  the  solid  phase  of  Che  system 
are  as  follows. 

1.  The  spheres  in  the  solid  phase  of  the  system  are  of  equal 
radii. 

2.  The  material  properties  of  the  solid  phase  are  constant. 

3.  The  spheres  in  the  solid  phase  are  arranged  in  ideal  packing 
configurations. 

With  these  restrictions  Eqn.  (4.20)  is  rewritten  as  follows. 


o  (1+v  )  C  b  ,  I  (b  ,)  +  3  C.  J 

s  s  [_  -2  k-1  -2  J  f  u 


~  (1-2V  )  Cs  b_2  +  Z  (b_2)k  +  Cfef 

s  k=l 


(4. 21) 


where  K  3  effective  bulk  modulus 

Cg  3  volume  fraction  of  the  solid  phase 
Cf  3  volume  fraction  of  the  fluid  phase 
5U  3  fluid  pressure 
ef  3  volumetric  strain  of  fluid 
Gs  3  shear  modulus  of  solid  phase 
-s  3  Poisson' 3  ratio  for  the  solid  phase 
M  3  the  number  of  contacts  on  one  sphere 
u_2  3  the  constant  of  superposition  resulting  from  a 

uniform  radial  pressure  on  the  surface  of  the  sphere 
(b_2)k  *  the  constant  of  superposition  for  contact  k. 


The  res cr i.c t ions  on  the  so L id- f Luid  system  requires  that  all  spheres 
are  in  static  equilibrium,  and  that  the  contacts  on  a  particular  sphere 
occur  in  pairs  along  an  axis  passing  through  the  origin  of  the  sphere. 
Furthermore,  all  spheres  contained  in  the  system  are  of  equal  radii  and 
have  the  same  material  properties.  These  restrictions  allow  the  use  of 
the  results  given  in  Section  3.4,  to  determine  the  constant  of  super¬ 
position,  The  results  contained  in  this  section  contain  values 

of  b_9  for  the  boundary  conditions  on  a  single  sphere  resulting  from 
Hertzian  contact  with  surrounding  spheres  and  an  uniform  radial  pres¬ 
sure  applied  on  the  surface  of  the  sphere.  The  effective  bulk  modulus, 

K,  given  by  Eqn.  (4.21)  will  be  de termined  .  for  the  cases  when  the 
displacements  or  stresses  are  specified  on  the  boundary  of  the 
representative  sample. 

4.2. 1.1  Displacements  Specified  on  the  Boundary  of  a  Representa¬ 
tive  Sample  of  a  Granular  System 

The  effective  bulk  modulus,  R,  will  be  determined  for  the  case  when  tne 
displacements  are  specified  on  the  boundary  of  a  representative  volume 
element.  The  rectangular  coordinate  system,  (x,  y,  z)  is  used  to  define  a 
point  in  the  representative  volume.  The  representative  volume  will  be 
suojecteo  to  a  macroscopicaiiy  homogeneous  displacement  field.  The 
displacement  fiexd  to  be  imposed  is  given  below 


u  (xj  =*  x  £ 

X  XX 


Jy  (y)  ’  y  £yy 


U  (z)  -  Z  e 

z  zz 


(4.22a) 


(4.22b) 


(4.22c) 


where 


=  the  displacement  in  the  x  coordinate  direction 


ux 

uv  =  the  displacement  in  the  y  coordinate  direction 
=  the  displacement  in  the  z  coordinate  direction 
x x  =  strain  in  the  x  coordinate  direction 
-  strain  in  the  y  coordinate  direction 
£2z  =  strain  in  the  z  coordinate  direction 

The  strain  field  appearing  in  Eqns.  (4.22)  is  assumed  to  be  constant 
with  respect  to  the  location  within  the  representative  volume.  The 
displacement  vector  at  any  point,  (x,  y,  z),  within  the  representative 
volume  is  given  by 
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where  u 


e 

e 

e 


x 

y 

z 


the  displacement  vector 

the  unit  vector  in  the  x  coordinate  direction 
the  unit  vector  in  the  y  coordinate  direction 
the  unit  vector  in  the  z  coordinate  direction 


Substitution  of  Eqns.  (4.22)  into  Eqn.  (4.23)  gives  the  following 
displacement  vector  in  terras  of  the  strain  field 


:xx  ex  +  y 


• yy 


czz  ez 


iT  =  x 


(4.24 


To  determine  tne  displacement  vector  at  a  point  of  contact  between  two 
spneres,  it  wiii  be  assumed  tnat  tms  vector  is  determined  by  evaluating  Eqn. 
(4.24)  at  tne  coordinates  of  tne  contact.  A  typical  sphere  contained  in  tne 
representati ve  volume  element,  witn  tne  origin  at  point  (x  .  yQ,  zQ) ,  is  snown 
m  Fig.  A. 2.  In  Fig.  4.2,  the  spnere  is  snown  in  contact  witn  two  adjacent 
spneres.  The  local  coordinate  directions,  x,  y,  z,  shown  in  Fig.  4.2  are 
parallel  to  the  x,  y  and  z  coordinate  directions,  respectively.  The  local 
coordinate  system  (x  ,  y  ,  z  ),  is  positioned  so  that  tne  z„  axis  is  directed 
througn  the  origin  of  the  spnere  and  contact  pair  m.  The  angles  6^  amd  ym 
define  the  position  of  the  local  coordinate  system  (x^,  ym,  zm) ,  relative  to 
tne  coordinate  system  (x,  y,  z).  The  contacts  occurring  along  the  zffl  axis  are 
labeled  as  contact  1  and  contact  2,  as  shown  in  Fig.  4.2.  The  coordinates  of 
tne  center  of  contact  1  in  reference  to  the  coordinate  system,  (x,  y,  z)  are 
given  by 


x  =  x  +  R  sm6  cosy 
1  o  m  m 


(4.25a) 


yi  “  yo  +  R  Sln*m  3in*m 


(4.25b) 


z.  *  z  +  R  cosy 
[  o  m 


(4.250 


wnere  x. 


-  x  coordinate 
»  y  coordinate 
*  z  coordinate 


of  tne  center 
of  tne  center 
of  tne  center 


of  contact  1 
of  contact  1 
of  contact  1 


R  -  radius  of  tne  spnere 


In  a  similar  manner,  Che  coordinates  of  Che  center  of  contact  2,  in 
reference  to  Che  coordinate  system,  ( x , y , z ) ,  are  given  by 


x2  =  xQ  -  R  sinS^  cos^ .  (4.26a) 

y2  =  y0  -  R  sinSm  sin^m .  (4.26b) 


z2  3  zo  "  R  cos  3m 

Using  Eqn.  (4.24),  the  displacement  vector  may  be  evaluated  for  the 
points  (x^,y^,z^)  and  '  The  displacement  vector  at  point  (x 

yi.Zf)  is  given  by 

U1  "  (xo  +  R  3inem  cos*m)  exx  ;x 
+  (yQ  +  R  sin6n  sin*ffl)  eyy  ey 

+  (z0  R  cos<|»m)  zzz  ez .  (4.27) 

The  displacement  vector  evaluated  at  point  ^x2’)r2,zz^  *3  8^ven  by 

u7  *  (x  -  R  sing  cos<i>  )  e  e 
4  o  m  m  xx  x 

>  (y0  -  R  sinBm  sin^)  eyy  ey 


>  (z0  -  R  cos^)  £zz  ez 


(4.28) 
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In  Eqns.  (4.27)  and  (4.28),  Che  subscripts  l  and  2  denote  Che  displace¬ 
ment  vectors  at  contacts  1  and  2,  respectively.  The  magnitudes  of  the 
displacements  in  the  positive  zm  coordinate  direction  may  be  deter¬ 
mined  for  contact,  1  and  2  by  taking  the  dot  product  of  the  displace¬ 
ment  vectors  at  these  contacts,  with  a  unit  vector  in  the  positive  zm 
coordiante  direction.  The  unit  vector  in  the  positive  zm  coordinate 
direction,  ezm,  is  given  by 

®zm  *  3inSra  cosK  +  sinSm  sinV  *  y  *  cos  Bra  *z  •  -(4-29) 

The  displacement  at  contact  1  in  the  positive  zra  coordinate  direction 
is  given  by 


“zm  *  “1  -  ®zm  *  Scx^o  +  R  sinBm  C0S^ra)sinBm  cos^m 

+  cyy(y0  +  R  3in6n  sin^)  sin&m  sini|/ffl 
*  ezz(zo  +  R  cosBrn)  cosBm  .  (4.30) 
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where  u^m  3  the  displacement  at  contact  1  in  the 


positive  zra  coordinate  direction. 


The  displacement  at  contact  2  in  the  positive  zm  coordinate  directior 
is  given  by 


exx  (*o  "  R  sin*m  cos^tn)  siniim  cos'V 


+  eyy  “  R  sinSm  sinV  sin3m  sinV 

+  £zz  (zo  “  R  c°*Sm>  cos3m .  (4.31) 

where  =*  the  displacement  at  contact  2  in  the  positive 

zm  coordinate  direction. 

Having  determined  the  displacements  in  the  positive  zm  coordinate 
directions  for  contacts  1  and  2,  the  change  in  the  diameter  of  the 
sphere  along  the  zm  coordinates  axis  is  given  by 

&uzra  3  uzm  "  uzm  *  2R  si-n23ra  cos2^m  +  2R  e  sin2S  sin2^ 

yy  m  i 

2R  ezz  cos2Sm .  (4.32) 

where  Auzm  *  t^ie  cbange  in  the  dimension  of  the  sphere  along 
the  axis  zra. 

The  Eqn.  (4.32)  gives  the  change  in  the  diameter  of  a  sphere  resulting 
from  contact  pair  ra,  directed  along  an  axis  passing  through  the  oiigin 
of  the  sphere,  as  shown  in  Fig.  4.2.  The  position  of  this  axis  with 
respect  to  a  global  coordinate  system,  (x,y,z),  is  defined  by  the 


angles  Sm  and  ^m.  If  Che  spheres  3hown  in  Fig.  4.2  are  all  of  Che 
same  radii  and  maCerial  propercies,  Chen  Che  deformaCion  occurring  aC 
Che  cenCer  of  each  contact  would  be  equal.  The  magnitude  of  this 
displacement  would  be  half  of  that  given  by  Eqn.  (4.32).  This  magni¬ 
tude  is  given  by 


=  R  (exx  sin2sm  cos2^m  + 


Eyy  sin23m  sin2*m  +  ez-  c°s23rj 


(4.33) 


where 


um  3  Che  displacement  at  the  center  of  either 
contact  for  contact  pair  m. 


The  direction  of  Che  displacement  given  in  Eqn.  (4.33)  is  toward  the 
origin  of  the  sphere  if  is  negative,  and  away  from  the  origin  of 
Che  sphere  if  u^  is  positive.  The  result  given  by  Eqn.  (4.33)  will 
be  used  to  determine  the  effective  bulk  modulus,  k,  for  the  case  when 

the  displacements  are  known  on  Che  boundaries  of  the  representative 
volume  element. 

For  the  case  when  the  displacements  are  known  on  the  boundaries  of 
the  representative  volume  element,  Eqn.  (4.21)  is  written  as 
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All  Che  Cerras  appearing  in  Eqn.  (4.34)  have  che  same  definitions  as  in 
Eqn.  (4.21),  except  that  the  superscript,  u,  has  been  used  on  the 
constants,  b_2  and  (b_2)m>  indicate  that  these  constants  are 
determined  for  the  case  when  surface  displacements  on  a  sphere  are 
specified.  From  Section  3.2,  the  constant  of  superposition,  b u ^  ,  is 
given  by 


( l-2v) 
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(4.35) 
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where  R  =  radius  of  sphere 

G  *  shear  modulus  of  sphere 
v  3  Poisson's  ratio  for  sphere. 

The  constant  appearing  in  Eqn.  (4.35)  is  determined  from  the 
displacement  boundary  conditions  on  a  sphere.  Two  eaaee-of  boundary 

conditions  on  a  sphere.  Two  cases  of  boundary  conditions  exist  on  the 

sphere.  The  first  is  the  displacements  resulting  from  contacts  with 

adjacent  spheres  and  the  second  due  to  the  displacements  from  a  uniform 

radial  fluid  pressure  acting  on  the  sphere. 

u 

The  constant  sQ  has  been  determined  for  the  displacement  boun¬ 
dary  conditions  resulting  from  the  contact  of  a  sphere  with  two  adja¬ 
cent  spheres,  along  an  axis  passing  through  the  origin  of  the  sphere. 
For  this  case  the  constant  i3 

,u  u  ;u 

5o  "  K  5o .  (4.36) 
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terras  K0,  and  are  given  by 


3(1-v2)R2F 
8  a3  E 


(4. 37) 


2  2 
-(2sin  b'-Dsin  <ja ' 


(l-cos% ' ) 


(4.38) 


where  F  =  the  force  transmitted  through  the  contact 

R  3  radius  of  sphere 
a  3  radius  of  the  contact  area 
4> '  3  sin-1  dr) 


The  Eqns.  (3.16a)  and  (3.20)  may  be  combined  to  determine  the  radius  of 
contact,  a.  For  spheres  of  equal  radii  and  material  properties,  the 
radius  of  contact,  a,  is  given  by 


3/  3(1-v2)RF 


(4.39) 


Substitution  of  Eqn.  (4.39)  into  Eqn.  (4,37)  yields  the  following  value 


for  ku 


(4.40) 


Combining  Eqns.  (4.35),  (4.36)  and  (4.40),  the  constant  of  superposi¬ 
tion,  b^ ,  can  be  determined  for  one  pair  of  contacts  on  a  sphere. 
This  result  is  given  by 


G  ^o 

b-2  "  "  2(l-2\>) 


(4.41) 


The  terra  ->Q  appearing  in  Eqn.  (4.41)  is  given  by  Eqn .  (4.38).  This 
terra  is  a  function  of  the  angle  defining  the  contact  area,  +>'.  The 
angle,  S',  may  be  related  to  the  displacement  at  the  center  of  the 
contact.  This  relationship  is  determined  from  Eqns.  (3.16)  and  is 
given  by 


=  sin’1  (f)  3  sin"1  ( J  |u|  /r)  ,  (0  <  S’  <  T/2)  .  .(4.42) 


where  u  =»  dispacement  at  the  center  of  the  contact  area 

R  3  the  radius  of  the  sphere. 


The  absolute  value  of  u  is  used  in  Eqn  (4.42)  since  >  0.  The  value  of  u  for 
a  pair  of  contacts  on  a  sphere  is  determined  from  Eqn.  (4.33).  This  value  of  u 

is  an  approximation  but  should  be  sufficient  for  the  determination  of 

the  effective  bulk  modulus,  k. 

u 

The  value  of  ^  for  the  surface  displacements  resulting  from  a 
uniformly  applied  radial  pressure  on  a  sphere  may  be  determined  from 
Section  3.4.3.  This  value  is  given  by 


(l-2v)R  a 
2G  ( l+'j ) 


(4.43) 


where  cu  =»  fluid  pressure  on  sphere 


Substitution  of  Eqn.  (4,43)  into  Eqn.  (4.35),  gives  the  constant  of 
.  .  u 

superposition,  b_2>  f°r  a  sphere  subjected  to  a  uniformly  applied 
radial  pressure.  This  constant  is  given  by 


Subs  t  i.  tu  c  Lon  of  Eqns.  (4.4L)  and  (4.44)  into  Eqn.  (4.34)  yields  Che 
following  for  Che  effecCive  bulk  modulus. 
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where 


k  3  effecCive  bulk  modulus 

C3  *  volume  fraction  of  solids  in  the  solid-fluid  system 
Cf  *  volume  fraction  of  fluid  in  the  solid-fluid  system 
Gg  *  shear  modulus  of  the  solid  phase 
vg  =  Poisson's  ratio  of  the  solid  phase 
au  3  fluid  pressure 

ef  *  volumetric  strain  of  the  fluid  phase. 

Eq  =  constant  determined  from  displacement  boundary 
conditions  on  sphere 

M  3  number  of  contacts  on  a  single  sphere. 


The  effective  bulk  modulus,  k,  given  in  Eqn.  (4.45)  is  for  use  when  Che 
displacements  are  specified  on  Che  boundaries  of  the  solid-fluid 
system.  The  simpliciations  used  to  arrive  at  Eqn.  (4.45)  from  Eqn. 
(4.20)  were  that  the  spheres  in  the  system  are  of  equal  radii  and  have 
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the  same  material  properties.  These  spheres  are  situated  in  ideal 
packing  configurations. 

To  evaluate  Eqn .  (4.45),  the  volume  fractions  and  material  proper- 

u 

ties  of  the  solid  and  fluid  phases  must  be  known.  The  term,  ,0,  must 
be  known  for  each  contact  pair  occurring  on  a  sphere  in  the  solid-fluid 
system.  This  constant  may  be  obtained  when  the 

displacement,  u,  of  the  center  of  the  contact  region  is  known.  To 
determine  the  displacement,  u,  Eqn.  (4.33)  is  used.  The  displacement, 
u,  is  the  displacement  normal  to  the  region  of  contact.  This  displace¬ 
ment  was  determined  by  treating  the  solid-fluid  system  as  a  continuum 
undergoing  macroscopically  homogeneous  displacements  of  the  form  given 
in  Eqn.  (4.42).  The  displacement  vector  at  a  point  in  this  continuum 
was  taken  to  be  that  occurring  on  a  contact  at  this  point.  The  magni¬ 
tude  of  the  displacement  vector  normal  to  the  region  of  contact  was 
taken  as  u.  This  value  of  u  is  an  approximation.  Any  displacements  on 
the  region  of  contact  which  are  tangential  to  this  region  have  been 
neg lec  ted . 

The  angles,  i  m  and  v'm,  defining  the  location  of  pairs  of 

contacts  with  respect  to  a  global  coordinate  system,  must  be  known  to 

r  u 

evaluate  the  sum  involving  the  term  -» 0- 

4.2.  T.  2  Surface  Tractions  Specified  on  the  Boundary  of  the 
Representative  Sample. 

The  effective  bulk  modulus,  k,  will  be  determined  for  the  case 
when  the  surface  tractions  are  specified  on  the  boundary  of  a  repre¬ 
sentative  volume  element.  The  representative  volume  element  will  be 
subjected  to  a  macroscopically  homogeneous  stress  field.  The  stress 
field  to  be  applied  to  the  representative  volume  is  given  by 


(  V  )  =  ** 

XX  'XX 


(4.46a) 
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where  ~Xx  =  Che  stress  normal  Co  Che  xz-plane 

3 yy  =  Che  stress  normal  Co  Che  xz-plane 

3 zz  =  Che  stress  normal  Co  Che  yz-plane 

The  stress  field  given  by  Eqns.  (4.46)  is  constanc  with  respect  Co 
position  within  the  representative  volume  element. 

An  arbibitrary  plane  cut  through  the  representative  volume  is 
shown  in  Fig.  4.3.  The  unit  normal  vector  to  this  plane  is  designated 
as  n.  The  angles  ax,ay  and  &z  define  the  position  of  the  unit  normal 
vector,  n,  to  Che  x,  y  and  z  axes,  respectively.  The  stress  vector,  ?, 
on  this  plane  is  given  by 


T  =  axx  cosCtx  ex  +  ayy  cos^  +  a zz  cosot z  ez  .  .  .  (4.47) 


where  ex  =  the  unit  vector  in  the  x  coordinate  direction 

Sy  -  Che  unit  vector  in  the  y  coordinate  direction 
ez  s  the  unit  normal  vector  in  the  z  coordinate  direction. 


The  unit  normal  veccor,  n,  to  Che  plane  shown  in  Fig.  4.3  is  given  by 


n  =  cos*x  ex  +  cosu  e„  +•  cos^,  e. 


(4.48) 


The  normal  stress  on  the  plane  is  determined  by  the  dot  product  or  tne 
stress  vector,  T,  with  the  unit  normal  vector,  n.  Performing  this 
operation,  the  normal  stress  is  given  by 

Jn  =  axx  cos2ax  ♦  ~yy  cos2ay  +  ~2Z  cos2^ . (4.4Q) 

The  normal  stress  given  by  Eqn.  (4.49)  was  determined  by  treating  the 
representative  volume  element  as  a  continuum  subjected  to  the  stress 

field  given  by  Eqn.  (4.46).  The  plane  intersecting  the  representative 
volume,  shown  in  Fig.  4.3,  will  contain  a  number  of  contact  regions 
present  in  the  solid-fluid  system.  These  contact  regions  are  formed  as 

the  spheres  in  this  system  are  compressed  together  under  the  action  of 

the  stress  field  given  by  Eqns.  (4.46).  Such  a  contact  region  is  shown 
in  Fig.  4.4.  The  normal  force  transmitted  through  this  contact  region 
will  be  determined  as  the  product  of  the  normal  stress  acting  on  the 
plane  containing  the  contact  region  and  an  area  contained  in-  the  plane. 

In  Fig.  4.4,  the  rectangular  coordinate  system,  ( x^, ,  ym,  zm)  ,  is  positioned 
such  th?1"  r.'ue  axis  is  directed  through  the  centers  of  the  contact 
regions  for  contact  pair  ra  and  the  origin  of  the  sphere.  The  coordi¬ 
nate  system  (x^ , y^ , zm)  is  defined  relative  to  the  coordinate  system 
(x,y,z)  by  the  angles  3m  and  vm.  The  coordinates,  x,y  and  z  are 
parallel  to  the  global  coordinates  x,  y  and  z,  respectively.  In  Fig. 

4.4,  the  plane  containing  one  of  the  contact  regions  for  contact  pair  m 

is  shown.  The  unit  normal  vector  to  this  plane  is  given  by 

n  *  sind  cosv_  e*  +  sin2m  sinv_  e  *  cosd_  e,.  .  .  .  (4.50) 
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where  ex  =  unit  vector  in  the  x  coordinate  direction 

ey  =  unit  vector  in  the  y  coordinate  direction 
e2  =  unit  vector  in  the  z  coordinate  direction. 


The  stress  vector  occurring  on  the  plane  shown  in  Fig.  4.4  is  given  by 


T  =  -  sine  cosv  e  +  3  sinC  sin>i  e  +  :  cose  e  .  .(4.51) 

xx  m  mx  vv  m  mv  zz  mz 


The  stress  components ,  ,  0yy  and  a^,  appearing  in  Eqn.  (4.51)  are 

chose  applied  to  che  boundaries  of  Che  representative  volume  element. 
The  normal  stress  on  the  plane  defined  by  the  unit  normal  vector 
appearing  in  Eqn.  v4.50)  is  determined  as  the  doc  product  of  the 
vectors  given  in  Eqns.  (4.50)  and  (4.51).  This  normal  stress  is  given 


•2,2  2  2 

=  3v„Sin;3  cos  j,  +  3  sin  3  sin  -p 

m  m  yy  m  i 


where 


+  o  cos  ^3  ....  (4.52) 

zz  m  J 


~  Che  normal  stress  on  the  plane  containing  Che 
contact  regions  for  contact  pair  m. 


The  rorce  transmitted  through  a  contact  region  contained  in  the  plane 

m  • 

sudjected  to  Che  normal  stress,  3  is  given  by 


r  m '  “  m  1  *■’  m  ^ 


=  ^n^m-  ^  A-n . . 


where 


area  of  Che  pLane  containing  the  contact  region 
which  transmits  normal  stress  through  the 
contact  region. 


m*  ****** **  * 


irur/*^  mi^J1  m  ihhwt> m  iwwiiTtfiTWiMM' 


162 


£ 

is 


I 


H5 


v. 

w; 


>: 


i 


fc. 

<*  , 

& 


/•j 

v; 


In  Eqn.  (4.53)  Fm  is  Che  force  cransmicced  through  a  contact  contained 

in  the  plane  defined  by  che  unit  normal  vector  given  in  Eqn.  (4.50). 

m 

The  normal  stress,  "n,  appearing  in  Eqn.  (4.53)  is  given  by  Eqn. 

(4.52).  To  determine  che  force,  Fm,  transmitted  through  a  pair  of 
contacts  ra,  it  will  be  assumed  chat  che  spheres  which  are  intersected 
by  a  particular  plane  carry  the  same  loads  as  the  contact  regions 
contained  in  che  place.  Therefore,  the  normal  stress  occurring  on  a 
plane  will  be  transmitted  evenly  through  the  spheres  and  contact 
regions  intersected  by  the  plane.  The  area  A^,  appearing  in  Eqn.  (4.53) 
is  given  by 
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From  Eqn.  (4.53),  the  force  transmitted  through  contact  A  is  given  by 
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The  normal  stress,  0^ ,  is  given  by  Eqn.  (4.52).  This  normal  stress 
depends  on  the  angles,  2m  and  -n  ,  defining  the  contact  orientation 
with  respect  to  the  global  coordinate  system,  (x,y,z).  A  particular 
sphere  in  the  system  under  consideration  will  normally  be  subjected  to 
many  contacts.  The  force,  Fra,  transmitted  through  a  contact  will 
depend  on  the  contact  orientation  and  the  area,  A^. 


The  value  of  Fm  is  an  approximation  to  the  actual  force  trans¬ 
mitted  through  a  contact.  The  contact  forces  occurring  in  many  systems 
of  spheres  situated  in  ideal  packing  configurations,  may  not  be  deter¬ 
mined  from  the  equations  of  static  equilibrium.  Due  to  this  the 
approximation  of  the  contact  force,  Fm,  given  by  Eqn .  (4.53)  will  be 
used.  This  approximation  will  be  used  to  determine  the  effective  bulk 
modulus,  R,  for  the  case  when  surface  tractions  are  known  on  the  repre- 
sentative  volume  element. 

For  the  case  when  surface  tractions  are  specified  on  the  boundar¬ 
ies  of  a  representative  volume  element,  Eqn.  (4.21)  is  written  as 

follows. 
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All  the  terras  appearing  in  Eqn.  (4.56)  have  the  same  meaning  as  those 
appearing  in  Eqn.  (4.21).  The  superscript,  o,  has  been  used  on  the 
constants,  b_  2  and  (b_2)^i  Co  indicate  that  these  constants  are  deter¬ 
mined  for  the  case  when  the  surface  tractions  have  been  specified  on 
the  boundaries  of  the  representative  volume  element.  In  Section  3.2, 
the  constant  of  superposition,  ba?,  was  determined.  This  constant  is 
given  by 
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(4.57) 


where 


v  ■ 


poison's  ratio  for  the  spheres. 


The  constant,  i  Q ,  appearing  in  Eqn.  (4.57)  is  determined  from  the 
surface  tractions  present  on  a  particular  sphere.  Two  types  of  boun¬ 
dary  conditions  will  exist  on  the  spheres  in  the  solid-fluid  system 
under  cons iderat ion.  The  first  of  these  boundary  conditions  are  the 
surface  tractions  resulting  from  contacts  with  adjacent  spheres.  The 
second  boundary  condition  is  due  to  a  uniform  fluid  pressure  acting  on 
the  surface  of  the  sphere. 

The  constant,  has  been  determined  for  the  surface  tractions 

which  result  when  a  sphere  is  in  contact  with  two  adjacent  spheres. 

The  contacts  made  with  the  two  adjacent  spheres  are  along  an  axis  which 
passes  through  the  center  of  the  sphere.  For  this  case  the  constant 
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The  terras  ku  and  appearing  in  Eqn. 


(4.45)  are  given  as  follows. 
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F  =  the  force  transmitted  through  the  contact 
R  *  radius  of  the  sphere 


a  =  radius  of  the  contact  area 


1C  5 


The  Eqn.  (4.39)  gives  Che  relationship  between  the  radius  of  the 
contact  region,  a,  and  the  force  transmitted  through  the  contact,  F. 
Using  this  relationship  the  constant,  ,  given  by  Eqn.  (4.59)  may  be 
rewritten  as 


2  E 


(4.61 ) 


n  (  1-v  ) 


where  E  =  modulus  of  elasticity  for  the  sphere 


Combining  Eqns.  (4.57),  (4.58),  and  (4.61),  the  constant  of  super- 

.  .  .  o  .  . 

position,  b_2>  1S  gi-ven  by 

b"  9  - . (4.62) 
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where  G  =  the  shear  modulus  of  the  sphere 

—  c 

The  quantity,  ?0,  appearing  in  Eqn.  (4.62)  is  given  by  Eqn.  (4.60). 

-a  . 

From  Eqn.  (4.60)  it  is  seen  that  is  a  function  of  the  angle,  i  , 
defining  the  contact  region.  This  angle  may  be  related  to  the  force, 

F,  transmitted  through  the  contact  through  Eqn.  (4.39).  This  relation¬ 
ship  is  g iven  by 
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The  approximation  given  by  Eqn.  (4.53)  will  be  used  to 

determine  F,  which  is  required  in  order  to  evaluate  the  effective  bulk 
modu  lus  ,  K. 

The  value  of  determined  for  the  case  when  a  uniform  rdial 

pressure  is  acting  on  the  surface  of  a  sphere  was  determined  in  Section 
3.4.3.  This  value  is  given  by 


’o  "  -u 


(4.64) 


where  =  fluid  pressure  on  the  sphere 


Substitution  of  Eqn.  (4.64)  into  Eqn.  (4.57)  gives  the  constant  of 
superposition,  b_7>  appearing  in  Eqn.  (4.56).  This  constant  is  given 
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(4.65) 


Substitution  of  Eqns.  (4.62)  and  (4.64)  into  Eqn.  (4.56)  yields  the 
following  expression  for  the  effective  bulk  modulus. 
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Where  K  =  effective  bulk  modulus 

C  =  volume  fraction  of  the  solids  in  the  solid-fluid  svstem 

3 

C.  =  volume  fraction  of  the  fluid  in  the  solid-fluid  system 
G  =  shear  modulus  of  the  solid  phase 
v  =  Poisson's  ratio  of  the  solid  phase 
=  fluid  pressure 

e.  =  volumetric  strain  of  the  fluid  phase 

=  constant  determined  from  the  stress  boundarv  conditions 

o 

on  a  single  sphere 

M  =  number  of  contacts  on  a  single  sphere. 

The  Eqn.  (4.65)  gives  the  effective  bulk  modulus,  K,  of  a  system 
comprised  of  a  number  of  equal  spheres,  in  ideal  packing  configurations, 
with  the  surrounding  void  space  being  filled  with  a  fluid  phase.  The 
effective  bulk  modulus,  K,  given  by  Eqn.  (4.65)  is  used  when  the  normal 
stresses  are  known  on  the  boundaries  of  a  representative  volume  element. 

To  use  Eqn.  (4.65),  the  volume  fractions  of  the  solid  and  liquid  phases 

present  in  the  system  must  be  known.  The  constants,  %  ,  must  be  deter- 

o 

mined  for  every  contact  pair  occurring  on  a  single  sphere  in  the  system. 

To  determine  these  constants  the  force,  F  ,  occurring  on  a  contact  pair 

m 

m  is  approximated  by  Eqn.  (4.53)  The  approximation  to  the  fore  e  >  ^m  >  was 

arrived  at  by  considering  Che  solid-fluid  system  to  be  a  continuum.  The  normal 

stress  acting  on  a  plane  containing  a  particular  contact  region  was 

determined  and  then  the  contact  was  said  to  carry  a  portion  of  this 

normal  stress.  Any  tangential  shear  stresses  acting  on  the  plane 

containing  the  contact  region  were  neglected  so  that  the  solutions 

given  in  Chapter  3  ,  could  be  used  m  determining  the  effective  bulk 

raodu lus . 
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4,2.2  Discussion  of  Che  Effective  3ulk  Modulus  Determined 
from  the  Volume  Averaged  Approach 

In  the  previous  two  subsections  of  this  chapter,  the  effective 
bulk  modulus,  K,  for  an  idealized  solid-fluid  system  was  determined 
using  a  volume  averaged  approach.  The  idealized  solid-liquid  system 
under  consideration  consists  of  a  number  of  spheres  in  contact  sur¬ 
rounded  by  a  fluid  phase.  In  order  to  obtain  a  simple  expression  for 
the  effective  bulk  modulus,  K,  of  such  a  system  some  restrictions  on 
the  geometry  and  material  properties  of  the  spheres  in  the  system  had 
to  be  imposed.  These  restrictions  are  as  follows: 

1.  All  the  spheres  in  the  system  have  the  same  material  proper¬ 
ties. 

2.  The  spheres  in  the  system  are  of  equal  radii. 

3.  The  spheres  in  the  system  exist  in  ideal  packing  configura¬ 
tions. 

The  first  of  these  restrictions  is  the  most  reasonable  with 
respect  to  modeling  actual  systems.  The  second  and  third  restrictions 
are  less  realistic  with  respect  to  modeling  actual  systems.  In  some 
cases  all  the  spheres  in  a  solid-fluid  system  may  be  of  equal  radii 
thus  making  the  second  restriction  valid.  The  third  restriction  which 
requires  the  spheres  in  the  solid-fluid  system  to  exist  in  ideal  pack¬ 
ing  configurations  is  the  most  unreasonable.  In  actual  systems  compos¬ 
ed  of  spheres  in  contact,  the  spheres  will  be  arranged  in  a  random 
order.  In  such  a  system  it  would  be  impossible  to  identify  the  loca¬ 
tions  of  the  contacts  present  on  a  single  sphere  so  that  statistical 
data  of  some  form  would  be  required  to  evaluate  the  effective  bulk 


modulus.  By  restricting  the  spheres  in  the  solid-fluid  system  to  be  of 


equal  radii  and  arranged  in  ideal  packing  configurations,  the  location 
of  the  contacts  on  a  single  sphere  witnin  the  svstem  may  be  identified. 
This  allows  for  the  determination  of  the  effective  bulk  modulus  without 
Che  use  of  statistical  data.  In  fact,  the  three  restrictions  thus  far 
mentioned  allows  the  determination  of  the  effective  bulk  modulus  by 
considering  only  one  sphere  and  ics  interactions  with  adjacent  spheres 
and  the  fluid  phase  present  in  the  solid-fluid  system. 

Other  restrictions  were  imposed  on  the  stresses  and  displacements 
occurring  on  the  regions  of  contact  between  adjacent  spheres.  These 
restrictions  are  as  follows. 

1.  The  contacts  on  a  particular  sphere  appear  in  pairs  with  each 
pair  occurring  along  an  axis  passing  through  the  origin  of  the 
s  phere. 

2.  The  stresses  and  displacements  occurring  on  a  particular 
region  of  contact  are  axisymme tr ic  with  respect  to  an  axis 
passing  through  the  center  and  perpendicular  to  the  region  of 
contac  t . 

3.  The  displacements  and  stresses  normal  to  the  region  of  contact 
may  be  determined  from  Hertz  contact  theory. 

The  purpose  of  these  restrictions  was  to  allow  the  use  of  the  solutions 
given  in  Chapter  3  ,  in  determining  the  effective  bulk  modulus.  Any 
non-symme tr ic  displacements  or  stresses,  tangential  to  the  region  of 
contact  have  been  neglected  in  the  determination  of  the  effective  bulk 
modulus  of  the  idealized  solid-fluid  system. 

Using  the  restrictions  thus  far  mentioned,  the  effective  bulk 
modulus,  K,  of  the  idealized  solid-fluid  system  was  determined  for  two 
types  ot  boundary  conditions  present  on  a  representative  volume  ele- 
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The  firsc  of  the  boundary  conditions  to  be  considered,  corre- 
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sponded  to  the  case  when  the  iisp  lac  emen  t  s  are  specified  on  the  surface 
of  the  representative  volume  element.  The  displacement  occurring  on  a 
particular  contact  was  approximated  by  treating  tne  solid-liquid  svstem 
as  a  continuum.  The  displacement  on  a  particular  contact  was  then 
determined  as  the  displacement  occurring  in  a  continuum  at  the  location 
of  the  contact.  This  type  of  approximation  is  used  because  the  actual 
value  of  the  displacement  on  a  particular  contact  is  indeterminate  for 
many  ideal  packing  configurations.  For  the  case  of  a  simple  cubic 
packing  configuration,  the  actual  displacements  occurring  on  the 
contacts  may  be  determined.  The  approximation  used  to  determine  the 
displacement  on  a  contact,  yields  the  actual  displacement  for  this 
case.  The  second  type  of  boundary  conditions  which  were  considered, 
corresponded  to  the  case  when  the  surface  tractions  are  specified  on 
the  surface  of  the  representative  volume  element.  For  this  case  the 
force  transmitted  through  a  particular  contact  was  again  approximated 
by  considering  the  solid-fluid  system  to  be  a  continuum.  The  normal 
stress  acting  over  a  portion  of  the  place  containing  a  particular 
contact  region  is  taken  to  be  the  force  transmitted  through  the 
contact.  This  approach  to  determine  Che  contact  forces  vields  the 
actual  value  of  this  torce  tor  the  case  of  a  simple  cubic  packing 
c  onf igurac ion . 

The  exoressions  for  the  effective  moduli  obtained  in  the  previous 
two  subsections  represent  very  idealized  conditions  as  apparent  from 
the  restrictions  and  approximations  used  in  Che  derivation.  In  using 
these  expressions  to  determine  the  effective  moduli  for  an  actual 
solid-fluid  system,  the  accuracy  of  the  value  obtained  for  the  effec- 
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tive  bulk  modulus  should  improve  as  this  system  becomes  better  repre¬ 


sented  by  the  idealized  solid-liquid  system.  Predictions  of  the  effec¬ 


tive  hulk  modulus  should  improve  for  hydrostatic  stress  states.  This 


is  because  there  would  be  no  shear  stresses  developed  in  the  idealized 


solid-liquid  system.  Any  shear  stresses  which  develop  due  to  non-hy¬ 


drostatic  stress  states  have  been  neglected  in  the  determination  of  the 


effective  moduli.  3v  the  same  reasoning,  the  predictions  of  the  effec¬ 


tive  bulk  modulus  should  improve  when  the  displacement  components  are 


all  of  equal  magnitude. 


4.3  Energy  Approach  for  Determining  Bounds  on 


the  Effective  Bulk  Modulus 


The  theorems  of  minimum  potential  energy  and  minimum  complementary 


energy  will  be  used  to  determine  upper  and  lower  bounds  on  the  effective 


bulk  modulus  of  an  idealized  solid-fluid  system.  These  theorems  and 


their  use  in  determining  upper  and  lower  bounds  for  the  effective  bulk 


modulus  will  be  discussed  below. 


To  discuss  the  theorems  of  minimum  potential  energy  and  minimum 


complementary  potential  energy,  an  elastic  body  in  static  equilibrium 


is  considered.  The  boundary  conditions  prescribed  on  such  a  body  are 


given  ov 


( 4 . 66a 


4 . 66b ) 


where  z.  .  *  stress  tensor 
4J 


u  =  components  of  displacement  vector 
i 

n.  =  components  of  the  unit  outward  normal  vector 

J 

on  the  surface 

I\  =  components  of  the  stress  vector  on  the  surface 
=  components  of  the  displacement  vector 
on  the  surface 

=  surface  over  which  T ^  is  prescribed 

S  =  surface  over  which  u.-  is  prescribed 
u  L 


Cartesian  tensor  notation  is  used  in  Eqns.  (4.66)  and  will  continue  to 
be  used  in  this  section.  The  potential  energy  functional  is  defined  as 
fo 1  Lows : 


a 


/ 


w  ( £ . . )  dV 

in 


dV 


ut  dS 


. (4.67) 


a 

where  *  components  of  the  body  force  vector. 

In  Eqn.  (4.67)  the  term  W  (e„)  is  called  the  strain  energy  function. 
This  function  is  given  by 


1/2  CijkZ  £ki£ij 


where  C..,  .  =  elastic  stiffness  tensor 
ijki 

,  ,  =  strain  tensor 

k^ 


(4.68) 


AdmissibLe  displacement  fields,  u^ ,  are  defined  as  any  continuous 
displacement  fields  satisfying  the  boundary  condition  given  ov  Eqn. 
(4.66b),  but  otherwise  arbitrary.  Under  these  restrictions  the  theorem 


■\. 


of  minimum  potential  energy  says  that  of  all  the  admissible  displace¬ 
ment  fields,  u^.  the  one  that  satisfies  the  equations  of  equilibrium 

» 

is  that  which  results  in  a  stationary  value  of  the  potential  energy 
functional.  This  may  be  written  as 


LL  =  (J,  >  0 


(4.69) 


where  U_  =  the  potential  energy  functional  evaluated  for 

any  admissible  displacement  field 
Ur  =  the  potential  energy  functional  evaluated  for  the 
displacement  field  which  satisfies  the  equations 


In  using  Eqn.  (4.69)  Co  determine  bounds  of  the  effective  bulk  modulus 
of  a  composite  rue  ter  ia  L ,  the  body  force  components,  F^,  will  be 
considered  to  be  zero  and  the  boundary  conditions  will  be  such  that 
displacements  are  prescribed  over  Che  entire  surface  of  the  body  under 
consideration.  For  this  special  case  the  potential  energy  functional 
given  by  Eqn.  (4.67)  reduces  to 


U 


W  (e.  )  dV 

4J 


(4. 70) 


The  Eqn.  (4.70)  gives  Che  strain  energy  of  the  body  under  considera¬ 
tion.  Therefore,  for  this  special  case  the  theorem  of  minimum  poten¬ 
tial  reduces  to  a  minimization  of  the  strain  energy  stored  in  the 


In  discussing  Che  theorem  of  minimum  complementary  energy,  an 
elastic  body  in  static  equilibrium  is  considered.  The  boundary  condi¬ 
tions  on  this  body  are  those  given  by  Eqns.  (4.66).  The  complementary 
.energy  functional  is  defined  as 


(4.71  ) 


The  term  W  (CT~„)  appearing  in  Eqn.  (4.71)  is  the  strain  energy  function 
expressed  in  terras  of  stresses.  This  function  is  given  by 


V  ( 


o  .  .) 

ij 


2  Sijk2.  akJl  °ij 


(4.72) 


where  S  =  L<=nsor  of  elastic  compliances 

The  stress  field,  CTij,  appearing  in  Eqn.  (4.71)  and  Eqn.  (4.72)  is 
admissible  when  the  boundary  conditions  given  by  Eqn.  (4.66a)  and  the 
equations  of  equilibrium,  are  satisfied.  Under  these  restrictions  the 
theorem  of  minimum  complementary  energy  says  that  of  all  the  admissible 
stress  fields,  the  one  which  satisfies  the  compatibility  equations  is 
that  which  results  in  a  stationary  value  of  the  complementary  energy 
functional.  This  statement  may  be  written  mathematically  as 


(3  -  U,.  >  0 

o  c  — 


(4.73) 


whe  re 


U „  =  the  complementary  energy  functional  evaluated 
for  any  admissible  stress  field 


r 


U,  =  che  complementary  energy  funccional  evaluated 
for  the  stress  field  which  satisfies  the 
compatibility  equations. 


In  using  Eqn.  (4.73)  to  determine  bounds  on  the  effective  bulk  modulus 
of  a  composite  material,  the  surface  tractions  will  be  prescribed  over 
the  entire  surface  of  the  body  under  consideration.  For  this  special 
case  the  complementary  energy  functional  is  given  by 


=  /  W  ( a  .  . )  dV 

a  1  ij 


(4. 74) 
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The  Eqn.  (4.74)  is  an  expression  for  the  strain  energy  of  the  body 
under  consideration.  Thus,  when  surface  tractions  are  prescribed  over 
the  entire  surface  of  the  body,  the  theorem  of  minimum  complementary 
energy  reduces  to  a  minimum  principle  for  the  strain  energy. 


li  V 

'v- 

« 

■?.  V, 
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4.3.1  Upper  and  Lower  Bounds  of  the 


Effective  Bulk  Modulus 


The  theorems  of  minimum  potential  energy  and  minimum  complementary 
energy  will  be  used  to  determine  upper  and  lower  bounds  on  the  effec¬ 
tive  bulk  modulus  for  a  composite  system.  The  Eqn.  (4.5)  may  be 


rewritten  as 


:  .  .  dV  =  —  <a .  . >  <c.  .>  V 
2  ij  lj 


(4.75) 
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where  <j ..>  =  volume  averaged  stress  tensor 

ij 

' ;  .  .  =  volume  averaged  strain  tensor 

ij 

■J .  .  =  infinitesimal  stress  tensor 
ij 

£  .  .  =  infinitesimal  strain  tensor 
i-J 

V  =  total  volume 


The  Eqn.  (4.75)  is  a  statement  that  the  strain  energy  stored  in  the 
different  phases  of  a  is  equal  to  the  strain  energy  expressed  in  terms 
of  volume  averaged  stress  and  strain.  The  integral  appearing  in  Eqn. 
(4.75)  is  taken  over  all  the  phases  contained  in  the  composite  mater¬ 
ial.  If  the  volume  averaged  stress  tensor,  <c ^ . > ,  or  the  volume  aver¬ 
aged  strain  tensor,  <£..>,  are  known  from  the  macroscopic  stress  or 

n 

strain  states  occurring  in  the  composite  material,  then  it  follows  from 
the  minimum  theorems  of  potential  and  complementary  energy  that 


1  (V  .£.  .  dV  >  <a.  .>  <e.  .>  V . (4.76) 

2  I  13  ij  ~  ij  ij 
V 

The  Eqn.  (4.76)  is  applicable  when  either  the  surface  displacements  or 
the  surface  tractions  have  been  specified  over  the  entire  boundary  of 
the  composite  material.  When  the  displacement  field  is  specified  on 
the  boundaries  of  the  composite  material,  then  an  admissible  strain 
field,  must  be  such  that  it  satisfies  the  displacement  boun¬ 

dary  conditions.  Then  the  stress  field,  J[j,  is  related  to  the 
strain  field,  —  ,  by 


rij  “  Cijk/£ii 


where  C 

ij  k 9. 


tensor  of  elastic  moduli. 


(4. 77  ) 
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When  the  surface  tractions  are  specified  on  the  boundaries  of  the 


romposite  material,  then  an  admissible  stress  field,  ;iji  must 


satisfy  these  stress  boundary  conditions.  Then  the  strain  field. 


is  related  to  the  infinitesimal  stress  field 


i  .  bv 


'  i  j  Si]kr'<:. 


(4. 78) 


where  S^.^,  =  tensor  of  elastic  compliances. 


An  upper  bound  on  the  effective  bulk  modulus,  K,  may  be  determined 


by  considering  the  composite  material  to  be  subjected  to  displacement 


boundary  conditions  which  result  in  a  uniform  volume  averaged  strain 


tensor,  <£„>.  The  displacement  boundary  conditions  of  this  type  are 


given  by 


<£kk>  - 
3  xi 


(4.79) 


where:  u^  =  surface  displacement  components 


<£kk> 


-  the  trace  of  the  volume  averaged  strain  tensor 


rectangular  Cartesian  coordinates  with 


respect  to  a  system  of  axes,  evaluated  on  the 


boundaries  of  the  composite  material 


The  volume  averaged  strain  field,  resulting  from  the  boun¬ 

dary  conditions  given  by  Eqn.  (4.79)  are  given  by 


1,4. 80  ) 


7TO 


where 


5. .  =  Kronecker  delta  = 
ij 


0,  i^j 


If  the  composite  material  is  an  isotropic,  linearly  elastic  material, 
then  the  volume  averaged  stress  field,  <j  — >,  occurring  in  the 
composite  material  is  given  by 


<0  .  >  =  K  <£..>  0-  . 
ij  ij  *-J 


(4.81  ) 


where 


K  =  the  effective  bulk  modulus 


Using  Eqns.  (4.76),  (4.80),  and  (4.81),  the  upper  bound  on  the  effec¬ 
tive  bulk  modulus  is  given  by 


K  <  Ku 


(4.82) 


where 


(4.83) 


To  evaluate  Eqn.  (4.83)  an  admissible  form  of  the  infinitesimal  strain 
tensor,  e^j,  is  required  for  all  phases  of  the  composite  material. 

The  infinites ima 1  stress  tensor,  o^j,  for  phase  k  is  given  by 


\  Ek 


(1+v.)  £iJ  +  ( 1+v  )  (l-2v  )  -KK  '  i- J 

K.  K.  k. 


(4.84) 


where  E^  -  the  elastic  modulus  of  phase  k 


'J >K  -  Poisson's  ratio  for  phase  k. 


A  lower  bound  on  the  effective  bulk  modulus,  K,  may  be  found  by 


t 


In  order  Co  evaluate  Eqn.  (4.88)  an  admissible  form  of  the  infinites¬ 
imal  stress  tensor  must  be  known  for  each  phase  contained  in  the  com¬ 
posite  material.  The  infinitesimal  strain  tensor,  e^j,  for  phase  k 
is  then  given  by 


(1+V  \ 

~E,  Jij  -  r  aKK  °ij  * 


(4.89) 


where  Ei,  =  the  elastic  modulus  of  phase  k 


v.  =  Poisson's  ratio  for  phase  k. 


To  summarize,  bounds  on  the  effective  bulk  modulus,  K,  are  given  by 


<  K  <  Ku 


(4.90) 


In  Eqn.  (4.90),  is  given  by  Eqn. (4. 88)  and  Ku  is  given  by  Eqn. 
(4.83). 


4.3.2  Upper  and  Lower  Bounds  on  the  Effective  Bulk 
Modulus  of  an  Idealized  Solid-Fluid  System 
The  general  form  of  the  upper  and  lower  bounds  on  the  effective 
bulk  modulus  have  been  determined  in  Section  4.3.1.  In  order  to  evalu¬ 
ate  the  expressions  for  the  upper  and  lower  bounds  the  phase  geometry 
and  phase  material  properties  must  be  specified.  The  remainder  of  this 
section  will  be  concerned  with  evaluating  the  expressions  for  the  upper 
and  lower  bounds  for  an  idealized  solid-fluid  system. 
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The  idealized  solid-fluid  system  to  he  considered  is  composed  of 
spheres  arranged  in  ideal  packing  configurations  and  surrounded  by  a 
fluid  phase.  The  spheres  in  this  assemble  will  all  have  the  same 
material  properties  and  will  be  of  equal  radii.  The  results  for  a 
sphere  subjected  to  axisymatric  surface  displacements  or  surface  trac¬ 
tions,  determined  in  Chapter  3  will  be  employed  to  determine  the 
bounds  on  the  effective  bulk  modulus.  Use  of  these  solutions  requires 
the  following  assumptions  concerning  the  contacts  between  adjacent 
spheres . 

1.  The  contacts  on  a  particular  sphere  occur  in  pairs  with  each 
contact  directed  along  an  axis  which  passes  through  the  origin 
of  the  sphere. 

2.  The  sphere  is  in  static  equilibrium  with  negligible  body 
forces . 

3.  No  non-axisymetr ic  surface  tractions  or  displacements,  occur 
on  the  region  of  contact  between  two  adjacent  spheres. 

The  first  of  these  assumptions  is  satisfied  automatically  by  consider¬ 
ing  only  spheres  in  ideal  packing  configurations.  The  geometry  for  a 
pair  of  contacts  is  shown  in  Fig.  4.4.  The  third  assumption  is  reason¬ 
able  since  the  bounds  on  the  effective  bulk  modulus  are  determined  for 
the  cases  of  the  representative  volume  element  undergoing  either  a  uni¬ 
form  strain  state  or  a  hydrostatic  stress  state.  For  the  case  of  the 
representative  volume  element  subjected  to  a  state  of  uniform  strain, 
Eqn.  (4.33)  gives  the  following  approximation  for  the  displacement,  u, 
occurring  at  the  center  of  a  contact  between  adjacent  spheres  contained 
in  the  system. 


(4.91  ) 


where  R  ~  Che  radius  Che  spheres  in  Che  solid-fluid  syscem 

<£,  >  =  Che  Crace  of  Che  volume  averaged  scrain  Censor. 


For  che  case  of  Che  represenCaC ive  volume  elemenc  subjecCed  Co  a  sCaCe 
of  hvdroscacic  scress,  Eqns .  (4.52  and  (4.53)  give  Che  following 
approximaC ion  for  Che  force  cransmiCCed  Chrough  a  conCacC  beCween  adja 
cenC  spheres  conCained  in  Che  syscem. 


<c  > 

F  =  — V .  (4-92) 


where  =  Che  area  which  Cransmics  Che  scress  normal  Co  Che 

conCacC  region  chrough  Che  conCacC. 

=  Che  Crace  of  Che  volume  averaged  scress  Censor. 


The  area,  Am,  appearing  in  Eqn .  (4.92)  will  depend  on  Che  packing 
con f iguraC ion  and  Che  oriencacion  of  Che  parCicular  conCacC  wich 
respecC  Co  Che  global  coordinaCe  sysCera  for  Che  represenCaC ive  volume 
elemenC.  To  simplify  Che  de Cermi na C i on  of  Che  bounds  on  Che  effecCive 
bulk  modulus  ic  will  be  assumed  ChaC  Che  force  acCing  on  all  conCacCs 
in  Che  sysCem  will  be  of  one  magnicude.  The  magnicude  of  Chis  force  i 
given  by 


<o  > 

p 


(4.93) 
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where  =  an  average  area  which  transmits  the  stress  normal 

to  the  contact  region  through  the  contact. 

The  approximations  to  the  displacement  at  the  center  of  a  contact 
and  the  force  transmitted  through  a  contact,  given  by  Eqn.  (4.91)  and 
Eqn.  (4.93),  respectively,  will  be  used  to  determine  the  bounds  on  the 
effective  bulk  modulus. 

To  determine  the  bounds  on  the  effective  bulk  modulus  the  integral 
appearing  in  Eqn.  (4.83)  and  Eqn.  (4.88)  must  be  determined.  The 
integral  appearing  in  Eqn.  (4.83)  is  used  to  determine  the  upper  bound, 
on  the  effective  bulk  modulus.  This  integral  must  be  determined  for 
the  case  when  the  displacements  are  specified  on  the  surface  of  the 
representative  volume  elements.  The  displacements  being  specified  on 
the  surface  of  the  representative  volume  element  are  given  by  Eqn. 
(4.79).  The  integral  appearing  in  Eqn.  (4.88)  is  used  to  determine  the 
lower  bound  on  the  effective  bulk  modulus.  This  integral  must  be 
determined  for  the  case  when  the  stresses  are  specified  on  the  surface 
of  the  representative  volume  element.  The  surface  tractions  to  be 
specified  are  given  by  Eqn.  (4.86).  The  integral  to  be  determined  can 
be  written  as 
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where  V  =  total  volume  of  Che  representative  sample 

Vs  =  volume  of  the  solid  phase 
Vf  =  volume  of  the  liquid  phase 
=  infinitesimal  stress  tensor 
£ i :  =  int ini tes ima 1  strain  tensor 
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The  superscripts,  s  and  f,  appearing  in  Eqn.  (4.94’>  are  used  to  denote 
quantities  for  the  solid  and  fluid  phases,  respectively.  From  the 
results  of  Section  4.3.1,  the  integral  taken  over  the  fluid  phase  is 
given  by 

u  ef  Vf . . (4.95  ) 


where  ee  3  the  volumetric  strain  of  the  fluid  phase 

-  =  the  fluid  pressure 


When  the  fluid  pressure,  ru,  is  known,  the  volumetric  strain  of  the 
fluid  phase,  ef,  is  determined  by 


The  tluid  pressure, 


appearing  in  Eqn .  (4.95)  and  Eqn.  (4.96)  is 


j!| 


liiri;ri;p 


taken  to  be  negative  for  gauge  pressures,  P,  greater  than  zero. 


The  solid  phase  of  the  representative  volume  element  will  be 


comprised  of  a  number  of  spheres  in  contact.  The  spheres  in  the  system 


bill  be  subjected  to  loads  resulting  from  contacts  with  adjacent 


spheres  as  well  as  from  the  fluid  pressure.  Using  the  principle  of 


superposition,  the  integral  taken  over  the  solid  phase  which  appears  in 


Eqn.  (4.97)  is  written  as 


. S  tS.  dV 


[a  .  .  +  o  .  F ] 
30  30 


+  ep. 


In  Eqn.  (4.97),  the  superscripts  c  and  p  are  used  to  denote  the  stress 


and  strain  fields  occurring  within  a  sphere,  which  result  from  contacts 


with  adjacent  spheres  and  the  fluid  pressure,  respectively.  The 


results  for  a  sphere  subjected  to  a  uniform  radial  pressure  are  given 


in  Section  3.4.3.  Using  this  solution,  and  -:P  •  are  given 


J  .  .  ~  J  y  ■  • 

1J  U  30 


(4.98) 


ij  ’  3  K  u  "ij 
s 


(4.99) 


where 


:u  =  the  fluid  pressure 


Ks  =  the  bulk  modulus  of  the  solid  phase 


O-  =  Kronecker  delta 


W 


Substitution  of  Eqns.  (4.98)  and  (4.99)  into  Eqn.  (4.94)  yields 
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(4.100) 
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The  terms  an ^  =4.  are  invariant  with  respect  to  coordinate  systems. 

Because  of  this  the  integrals  containing  these  terms  may  be  performed 
without  transforming  the  solutions  for  the  stress  and  strain  fields  due 
to  a  particular  contact  into  one  global  coordinate  system.  The  inte¬ 
gral  containing  these  terms  were  determined  in  Section  4.2.1.  The 
values  of  these  integrals  are  given  by 


f  dV  =  3  (1+v  )  C  M  b  _  V . (4.101) 
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where  M  =  the  number  of  contacts  between  a  particular 

sphere  and  adjacent  spheres 
Cs  =  the  volume  fraction  of  the  solids  contained 
in  the  system 

G„  =  the  shear  modulus  of  the  solid  phase 


Following  the  above  approach  would  yield  for  energy  functionals 


Uc  (M'  V 
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(4.103) 


where  M  -  the  number  of  contacts  on  a  particular  sphere 

“  the  an8les  defining  the  position  of  contact  m  with  respect 
to  a  global  coordinate  system 
9*  =  the  angle  defining  the  region  of  contact  on  a  sphere 
v  =  Poisson's  ratio  for  the  sphere 


For  spheres  arranged  in  a  particular  ideal  packing  configuration,  the  number 
of  contacts  on  a  particular  sphere,  N,  and  the  angles  and  9m,  are  known. 
Therefore  for  a  particular  ideal  packing  configuration,  the  quantities,  UU  and 
U°,  may  be  presented  as  a  family  of  curves  for  different  values  of  <p’  and  v  . 

w  ‘3 

The  expressions  for  the  upper  and  lower  bounds  on  the  effective  bulk 
modulus,  K,  may  now  be  determined.  The  upper  bound  on  the  effective  bulk 
modulus,  R^,  is  determined  for  the  case  when  the  displacements  are  specified 
on  the  representative  volume  element.  Combining  Eqns.  (4.41),  (4.83).  4.94), 
4.95),  (4.100),  (4.101),  and  (4.102) 
gives  the 

-  K1U 
K  -  tt- 
u  K2u 

The  values  of  and  K 2  are  given  by 


(4.104) 
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KlU=  4(  l+vs )  Cg  au  3  Ks  M  £  ♦  2  au 
2  Kgef  +  8(l+vg)  au  ef 
+  3(1-2VS)  cs  Kg  ucU  (4.105) 


K2  =  8(lk+Vg)  Cg  <£kk>2 


(4.106) 


The  lower  bound  on  the  effective  bulk  modulus,  K^,  is  determined  from  the  case 
when  the  surface  tractions  are  specified  on  the  representative  volume. 
Combining  Eqs.  (4.62),  (4.88),  (4.94),  (4.95),  (4.100),  (4.101),  and  (4.102) 
gives  the  following  for  , 


K2L 


r»L  rjL 

The  values  of  and  are  given  by 
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♦  24  (l-2vg)  (l+Vg)  Cg  Kg'  u/ 


In  summary,  the  upper  and  lower  bounds  of  the  effective  bulk  modulus,  ft, 
of  the  idealized  solidnfluid  system  under  consideration  are  given  by 

ft.  -  £  ft  S  R  (4.110) 

Li  li 

The  quantities,  ft.  and  ft  ,  are  determined  from  Eqns.  (4.101  to  4.106).  The 
expressions  for  the  upper  and  lower  bounds  show  that  the  voids  in  the  system 
have  an  effect  on  the  effective  bulk  modulus  when  the  fluid  pressure  is 
atmospheric.  This  was  not  the  case  for  the  effective  bulk  modulus  determined 
from  the  volumetric  averaging  approach. 


4.4  Effective  Poisson’s  Ratio 

In  this  section  an  effective  Poisson's  ratio,  v,  will  be  determined  for 
the  idealized  solidnfluid  system.  This  effective  Poisson's  ratio  will  be 
determined  through  the  definition  of  this  quantity.  Poisson's  ratio  is 
defined  as  the  negative  of  the  ratio  of  lateral  strain  to  longitudinal  strain 
for  a  sample  loaded  along  the  longitudinal  axis.  Fig.  4.5  shows  a  sample 
undergoing  uniaxial  loading  along  the  z  axis.  The  dotted  line  shows  the 
deformed  shape  of  the  sample  due  to  the  applied  load.  For  the  sample  shown  in 
Fig.  4.5,  Poisson's  ratio  is  determined  as 


v 


A  /L 

jl  y 

A  /L 
X  X 


(4.111) 


where  v  -  Poisson's  ratio 


The  effective  Poisson's  ratio  for  the  solid-’fluid  system  can  be  determined  in 


the  same  manner  as  that  given  by  Eqn.  (4.111).  The  idealized  solid-fluid 
system  consists  of  a  number  of  equal  spheres  arranged  in  ideal  packing 
configurations.  These  spheres  are  surrounded  by  a  fluid  phase.  The  following 
assumptions  are  required  in  order  to  allow  the  determination  of  the  effective 
Poisson’s  ratio, 

1 .  The  spheres  in  the  system  do  not  experience  any  rigid  body  motion. 

2.  Only  normal  forces  are  transmitted  through  contacts  between  adjacent 
spheres. 

3.  The  fluid  pressure  in  the  system  is  at  steady  state  conditions. 

The  first  assumption  is  necessary  so  that  the  geometry  of  the  spheres  in  the 

system  remains  known.  The  second  assumption  allows  the  use  of  the  results 

given  in  Chapter  3  in  determining  the  effective  Poisson's  ratio.  The  reason 

for  the  last  assumption  will  be  explained  presently. 

To  determine  the  effective  Poisson's  ratio,  the  idealized  solid-;fluid 

system  will  be  loaded  in  one  direction  only  as  shown  in  Fig.  4.6.  The  system 

may  initially  be  stressed  due  to  body  forces.  The  effective  Poisson's  ratio 

will  be  determined  from  the  displacements  occurring  in  the  system  when  the 

system  is  disturbed  from  its  initial  state.  The  addition  of  the  stress,  Aa  , 

zz 

as  shown  in  Fig.  4.6  may  cause  an  increase  in  the  fluid  pressure  as  well  as 

the  stresses  present  in  the  spheres  making  up  the  solid  phase.  The  effective 

Poisson's  ratio  of  the  system  may  be  determined  for  the  case  when  the  fluid 

pressure  has  dissipated  to  its  initial  value  (steady  state  condition).  For 

this  condition  the  additional  stress,  Ao  ,  will  be  transferred  to  the  solid 

zz 

phase  and  the  effective  Poisson's  ratio  may  be  determined  by  only  considering 
this  phase.  If  non-steady  state  conditions  were  considered,  a  sample  of  the 
material  from  within  the  system  would  be  subjected  to  multinaxial  stresses  in 
reference  to  its  initial  state.  Thus,  the  definition  of  Poisson's  ratio  could 
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not  be  used  to  determine  the  effective  Poisson's  ratio.  For  steady  state 
conditions,  the  change  in  stress  on  the  system  will  be  one  directional  and  the 
definition  of  Poisson's  ratio  may  be  used  to  determine  the  effective  Poisson's 
ratio.  Under  the  loading  shown  in  Fig.  4.6,  the  two  planes,  A-A  and  A'-.A', 
will  approach  each  other  by  some  amount,  while  the  two  planes,  BnB  and  B'-B', 
will  move  away  from  each  other  by  some  amount.  These  displacements  may  be 
determined  from  the  results  given  in  Chapter  3-  From  these  displacements  the 
effective  Poisson's  ratio  may  be  determined.  The  change  in  the  displacement 
vector  at  any  point,  (x,  y,  z),  in  the  sphere,  due  to  the  addition  load  is 
given  by 


rui  * 


z  T 

E  (TrJ_  ( 
m=l 
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-  iuli  J 


(x,y,z) 


(4.112) 


where  M  »  the  number  of  contacts  transmitting  force 
m  -  the  contact  of  interest 


In  Eqn.  (4.112),  the  superscripts  F  and  I  denote  final  and  initial  values  of 
the  displacement  vector,  in  reference  to  a  local  spherical  coordinate  system, 
respectively.  The  symbol,  A,  appearing  in  this  equation  is  used  to  denote  the 
change  in  the  displacement  vector.  The  transformation  matrices,  [T„]  and 

U 

[T,],  are  given  by  Eqns.  (3.140)  and  (3.141),  respectively.  The  displacement 

U  'S 

vector,  {Au} ,  in  reference  to  the  global  coordinate  system  is  given  by 


1  Au  1 


Au, 

Au, 


(4.113) 
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The  Eqn.  (4.112)  may  be  evaluated  at  contacts  contained  on  the  planes,  AnA  and 
®  A'-A'f  to  determine  the  strain  occurring  between  these  planes.  This  strain 

$  is  given  by 

(4.114) 
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a  contact  which  transmits  force. 

an  angle  defining  the  position  of  a  contact  transmitting 
force  to  the  efVtfce' 
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CHAPTER  5 

EFFECTIVE  MODULI  OF  AN  IDEALIZED  3~PHASE  SYSTEM 


5.1  Introduction 

In  the  previous  chapter  effective  moduli  for  an  idealized  soil-water 
system  was  derived.  This  model  may  be  applicable  for  a  fully  saturated 
system.  In  a  partially  saturated  system  where  soil,  water  and  air  are  in 
contact,  previously  derived  moduli  cannot  be  used.  The  compressible  nature  of 
the  gaseous  phase  plays  an  important  role.  In  this  chapter  an  attempt  is  made 
to  arrive  at  an  effective  moduli  for  the  soil-swatemair  3-phase  system  under 
idealized  conditions.  The  approach  used  by  Chang  &  Duncan  (1974)  is  used  in 
arriving  at  a  compressibility  parameter  for  this  air-swater  system. 

In  this  chapter,  in  the  development  of  the  effective  moduli  for  a  Biphase 
system  it  is  also  assumed  that  the  effective  stress  variation  in  a  partially 
saturated  granular  void,  changes  with  the  degree  of  saturation. 

5.2  The  compressibility  of  the  homogenized  pore  fluid  of  partially  saturated 

soil 

The  compressibility  of  the  horaogenzied  pore  fluid  1/Kp  is  defined  as  the 
volumetric  strain  induced  by  a  unit  change  in  pore  fluid  pressure.  It  can- be 
defined  by  the  following  equation 

,  d  e 

!_  _ p 

K  "  ~  d  o  5,1 

P  P 

where  e  -  volumetric  strain  of  pore  fluid 
P 
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0^  <*  pore  fluid  pressure 


The  homogenized  pore  fluid  pressure  can  be  expressed  as  a  function  of  air 
pressure,  water  pressure  and  surface  tension.  It  may  be  expressed 


o  =  f(o  ,o  ,  T  ) 
p  awe 


where  =  air  pressure  in  the  pore 

o^  *  water  pressure 

T  »  surface  tension 

c 


The  compressibility  of  tne  homogenized  pore  fluid,  therefore,  can  be  evaluated 
using  the  following  eq. 


-  ,*J-  3Tc 

3  a  3e  +  3  o  3e  +  3T  3  e  ^ 
a  p  w  p  c  p 


If  the  form  of  Equation  5.3  is  known,  the  compressibility  of  the 
homogenized  pore  fluid  can  be  evaluated  by  knowing  the  changes  in  air 
pressure,  water  pressure,  and  tension  due  to  unit  change  in  volumetric  strain. 
For  the  case  in  which  the  air  bubbles  are  occluded,  Eq.  5.3  reduced  to: 


For  this  important  practical  case,  the  compressibility  of  the  homogenized  pore 
fluid  can  be  evaluated  using  Boyle's  law  and  Henry's  law.  This  darivation  is 
given  in  the  following  section. 


Note  that  for  clays  with  low  degrees  of  saturation,  open  channels  are  likely 
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to  be  present.  The  behavior  in  such  a  case  is  extremely  complicated.  The 
approach  here  considers  only  granular  soils. 


5 . 3  Derivation  of  compressibility  expression: 

This  derivation  is  based  on  Boyle's  law  and  Henry's  law. 


Boyle's  Law 

Boyle's  Law  states  that  the  product  of  pressure  and  volume  of  a  gas  is 
constant  under  constant  temperature  conditions: 


Henry's  Law 

Henry's  Law  states  that  at  a  constant  temperature,  the  weight  of  gas 
which  can  be  dissolved  in  a  given  volume  of  liquid  is  directly 
proportional  to  the  gas  pressure. 


Let  Vd  be  the  volume  that  would  be  occupied  by  the  dissolved  gas  if 
it  was  extracted  from  the  liquid  and  compressed  to  the  pressure  acting  on  the 
fluid.  According  to  Boyle's  and  Henry's  Laws  this  volume  will  be  constant  and 
tnus  independent  of  the  pressure.  is  only  dependent  on  the  volume  of  water 
and  can  be  calculated  as  follows: 


V  ■  H  V 
d  w 


where  Vw  is  the  volume  of  water,  H  is  the  coefficient  of  solubility.  At  20°C, 
H  -  0.02. 


According  to  Dalton's  division  law,  the  saturated  water  vapor  pressure  in 
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che  free  air  does  not  obey  Boyle's  Law.  However,  the  saturated  water  vapor 
pressure  is  usually  very  small  (Schuurman,  1966)  and  tne  influence  can  be 
disregarded. 


According  to  Henry's  Law,  the  volume  of  air  dissolved  in  the  water  is 
proportional  to  the  volume  of  water.  The  rate  at  which  the  air  dissolves  in 
the  water  depends  on  the  air  pressure.  The  time  t  necessary  to  dissolve  the 
air  for  a  unit  change  in  pressure  was  shown  by  Beek  (1963)  to  be 


2 

r 


5.6 


where  Df  is  the  diffusion  coefficient  and  r  is  the  radius  of  the  bubble.  At 

^5 _ 2 


20°C,  Df  is  equal  to  10  cm  /sec.  For  small  air  bubbles,  the  time  necessary 
to  dissolve  the  air  is  very  small.  Therefore,  it  can  be  assumed  that  the  air 
dissolves  in  the  water  instantaneously  and  that  the  water  is  always  saturated 
with  dissolved  air. 


As  a  result  of  the  surface  tension,  T  ,  the  air  pressure  and  the  water 


pressure  of  an  airnwater  mixture  are  not  the  same.  The  value  of  o  and  a  are 

w  a 


related  by 


o  n  o 

a  c 


5.7 


where  a q  is  the  capillary  pressure  due  to  surface  tension.  oc  is  usually 


expressed  in  terms  of  the  meniscus  radii,  r  ,  r^,  uj  shown  in  Fig.  5.1. 
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In  a  mixture  of  air  and  water,  the  bubbles  do  not  necessarily  all  have  the 
same  shape  and  size.  Therefore,  it  is  assumed  that: 

T 

a  =*  2  5.9 

c  R 

c 

where  Rc  is  the  average  capillary  radius  for  the  mixture  of  air  and  water. 
Generally,  for  soils  with  high  degrees  of  saturation,  the  effect  of  oc  is 
small. 

If  e  is  the  volume  of  free  air  and  e  is  the  volume  of  dissolved  air  in 
a  s 

the  air~water  mixture,  then  it  follows  from  Boyle's  and  Henry's  Law  that 


lea  *  e3>  (pa  * 


°a> 


(eao  *  V  <Pa  * 


0  ) 
ao 


5.10 


where  p^  is  atmospheric  pressure,  eaQ  is  the  initial  air  volume  and  oaQ  is  the 
initial  air  gauge  pressure.  Equation  5.10  can  also  be  expressed  as: 


(e  h  e  )P  +■  (e„„  +  e  )a 
ao  a  a  ao  a  ao 

~  e  7~e 
a  s 


5.11 


Using  expressions  5.7  and  5.11,  the  pore  water  pressure  may  be  expressed  as: 


(e  n  e  )p  +  (e  „  e 
ao  a  a  ao  - 


s)  ao 


2T 


e  +  e 
a  s 


5.12 


The  compressibility  of  the  airnwater  mixture  is, 


202 


1 

iT 


de 


,  d 
1  ea 
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where  e  is  the  initial  voids  ratio, 
o 


do  (e  +  e  )P  +  (e  +  e  )o  2T  de 

_ w  ao  s  a _ ao  s  ao  _ c  _ c 
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5.13 


5.14 


In  terras  of  degree  of  saturation  S  and  void  ratio  e 


do 

de 


(e  -.*s  e  +HS  e  )  (u  +  P  ) 
o  oo  oo  ao  a 

(e  -  Se  +  HSe)2 

2T  dR 
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5.15 


Substituting  Eq.  5.13  into  Eq.  5.15,  the  compressibility  of  the  homogenized 
pore  fluid  can  be  expressed  as: 

1 


1_ 

K 


( 1 +eo) 


5.16 


(e  -*s  e  +HS  e  (o  +P  )  n  2t  dR  _ 
o  oo  ooaoa  c  cS 


(e-Se+HSe)' 
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2  dS  e 


where 


o 

e 

S 

H 


ao 


Bulk  modulus  of  the  airnwater  moisture 

initial  voids  ratio 

initial  degree  of  saturation 

voids  rat*o 

degree  of  saturation 

Coeff.  of  solubility 

Surface  tension 

Average  capillary  radius 

initial  gage  air  pressure 
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If  the  surface  tension  term  in  Eq.  16  is  neglected,  the  equation  reduces 
to  a  form  similar  to  that  suggested  by  Hilf  (1948)  and  Skempton  and  Bishop 
(1954).  The  comparison  between  test  data  obtained  by  Mitchell  et  al.,  (1965) 
and  the  results  predicted  by  Hilf's  and  Skempton  and  Bishop's  equation  is 
shown  in  Fig.  5.2. 

Schuurman  (1966)  took  surface  tension  into  account  by  assuming  that  the 
air  bubbles  were  all  of  equal  size  and  spherical  in  shape.  By  this 
assumption,  the  value  of  Rc  in  Eq.  5.9  would  be  equal  to  the  radius  of  the  air 
bubbles.  Taking  the  initial  radius  of  air  bubbles  as  RQ,  the  radius  after  a 
cnange  in  void  ratio  can  be  expressed  as 

e  1/3 
c  °  eaoJ 

By  substituting  this  expression  into  Eq.  5.15,  the  compressibility  of  this 
air-iwater  mixture  can  be  calculated.  A  comparison  between  test  data  obtained 
by  Mitchell  et  al.  (1965)  and  the  results  predicted  by  Schuurman's  equation  is 
shown  in  Fig.  5.2. 

It  seems  likely  that  the  discrepancies  between  experimental  data  and  the 
results  predicted  by  Hilf’s  and  Skempton  and  Bishop's  equation  are  mainly  due 
to  the  negligenece  of  the  effect  of  surface  tension.  The  discrepancies 
between  the  experimental  data  and  the  results  predicted  by  Schuurman's 
equation  are  probably  due  to  the  assumption  that  all  the  bubbles  are  spherical 
and  uniform,  and  that  the  total  number  of  bubbles  does  not  change  during  the 
compression  of  the  soil. 
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A  further  consideration  relative  to  surface  tension  effects  is  that  the 
capillary  pressure  does  not  necessarily  arise  only  from  the  bubbles.  It  can 
also  arise  from  the  other  water  surfaces  exposed  to  air  at  the  specimen 
boundaries.  Therefore,  it  seems  more  reasonable  to  use  an  empirically 
determined  average  capillary  radius  Rq  to  simulate  capillary  effects. 

It  has  been  found  that  the  calculated  results  agree  quite  well  with  the 
test  by  Mitchell  et  al.  (1965)  if  it  is  assumed  that  the  average  capillary 
radius  Rq  and  the  degree  of  saturation  are  related  as  follows: 


wnere  R  is  the  capillary  radius  at  saturation.  The  comparison  is  shown  in 
cs 

fig.  5.3.  Sf  is  the  lowest  degree  of  saturation  at  which  the  water  begins  to 

flow  freely.  Sf  may  be  found  from  permeability  tests  for  samples  with 

different  degrees  of  saturation.  For  the  clay  tested  by  Mitchell,  et  al.,  S^. 

was  found  to  be  zero.  It  was  also  found  that  the  relation  between  the  water 

pressure  and  the  degree  of  saturation  was  approximated  best  by  adopting  a 

value  of  2T  /R  »  5  psi.  At  2Q°C,  the  value  of  T  is  74  x  10^  kg/cra.  The 
c  cs  c 

value  of  2Tc/Rg3  -  5  psi  corresponds  to  an  effective  value  of  Rcg  equal  to 
4.23  x  10“4  cm. 

5 .4  Effective  Bulk  modulus  of  a  3~pbase  system. 

Previously  an  expression  was  derived  for  the  effective  bulk  modulus  of  a 
saturated  2-phase  soil-water  system.  The  bulk  modulus  expression  in  Eqn.  4.66 
is  further  modified  here  for  a  3-:phase  soilnwater-air  system  in  light  of  the 


(Data  from  Mitchell,  et  al.,  1965 
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compressibility  expression  for  a  waternair  mixture  derived  before. 


In  Eq.  4.66  the  parameters  affected  by  the  degree  of  saturation  are  c^., 

V  ef  and  Co* 


Assuming  that  the  pore  pressure  in  partially  saturated  systems  is  directly 
related  to  the  degree  of  saturation  via  X  .where  X  is  an  empirical  parameter 
given  in  Fig.  5.4  in  the  relation 


a*  =  a  -  cTa  +  x  (<*a  -  *,) 


The  volumetric  strain  for  a  partially  saturated  system,  ef,  can  be  determined 
from  the  following  expression 


'f  IS. 


where  K  -  Bulk  modulus  of  the  air-water  mixture. 
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CHAPTER  6 

CONCLUSIONS  AND  RECOMMENDATIONS 


From  this  investigation,  the  following  conclusions  can  be  drawn: 

1 .  The  modelling  of  even  saturated  systems  using  the  micromechanics  approaches 
is  quite  complicated,  and  the  degree  of  complexity  is  far  greater  with  such 
approaches  for  partly  saturated  systems. 

2.  Numerous  assumptions  need  to  be  made  in  order  to  arrive  at  the  stressn 
strain  response  of  partly  saturated  systems  using  the  raicromechanicsnbased 
models.  However,  one  could  obtain  some  insight  into  the  fundamental  behavior 
of  partly  saturated  systems  in  this  methodology. 

3.  The  expressions  developed  for  effective  moduli,  and  the  procedure  outlined 
for  poisson's  ratio  should  be  valuable  in  modelling  behavior  of  partly 
saturated  systems. 

4.  Although  the  immediate  usefulness  of  this  type  of  modelling  approach  may 
not  be  that  evident,  the  long  term  potential  for  improving  the  understanding 
of  the  behavior  of  partly  saturated  soils  is  great. 

The  recommendations  for  futher  research  include  the  possible  elimination  of 
some  of  the  assumptions  based  on  experimental  data,  and  validation  of 
predictions  by  the  newly  developed  models  on  various  partly  saturated  systems 
using  results  from  controlled  experiments.  Another  useful  further  research 
effort  would  be  to  extend  this  work  for  timendependent  behavior  and  strain^ 
rate  dependent  behavior. 
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